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1 Introduction 


The simulation of scattering events according to field theoretical models is at the core of 
particle physics, as it allows for both discoveries of particles and the validation of theories. 
The source of practical difficulties lies primarily in the presence of strong interactions. 
Indeed, if the initial and/or final states involve hadrons, then non-perturbative effects will 
have to be included at some point. It is only due to factorization that it is possible to 
obtain an accurate description using perturbative methods. For some observables, the 
evolution of the state of the system between the high scale, where fixed order perturbation 
theory is adequate, and the non-perturbative scale, where it is not, requires resummation 
of logarithmic enhancements due to soft and collinear singularities. This can be achieved 
either by analytic resummation methods, or by means of parton showers. The latter are 
the versatile tool of choice in the communication between theory and experiment. 

Parton showers have a long history [1-6]. They are based on a semi-classical picture 
of partons splitting into pairs of other partons. This picture corresponds to the collinear 
limit of QCD amplitudes in next-to-leading order (NLO) configurations. Thus, the lead¬ 
ing singularity of collinear-soft origin may be correctly reproduced. In consequence, the 
approach has leading logarithmic accuracy in the most general case of many partons in 
the high scale process. Including pure soft effects is necessary to correctly reproduce the 
first subleading logarithms. Recent years have seen quite some activity in this direction [7- 
13]. In particular, it has been possible to include leading colour soft effects using parton 
showers constructed around known subtraction schemes, such as the one of Catani and 
Seymour [14, 15]. In a series of papers [16-25], Nagy and Soper have proposed a different 
concept for a parton shower. Their construction should be able to include soft effects at 
subleading colour. Recently, they have also provided a public code [24], DEDUCTOR, that 
implements a large part of their idea, albeit with the exlusion of exact colour treatment, 
where an extension of the leading colour approximation is used, and spin correlations. 

The purpose of the present publication is to specify a matching procedure between the 
Nagy-Soper parton shower and a fixed order calculation at NLO. This can be achieved using 
at least two different ideas: Mc@Nlo [26, 27] and POWHEG [28-30]. The first method, 
which is the one we have chosen, consists in removing double counting contributions by 
expanding the parton shower to first order in the strong coupling and compensating for the 
terms, which are already present at fixed order. This is closely related to the possibility of 
obtaining a subtraction scheme from the parton shower, which was exploited in Refs. [31- 
34] for the case of the Nagy-Soper construction. In particular, the implementation within 
the HELAC-DIPOLES package [35] has proven invaluable as the basis for the present work. 

Having a working tool to produce parton shower matched event samples with inclusive 
next-to-leading order accuracy is necessary in order to assess the differences of the new 
shower concept and more established alternatives. This is a second goal of our publication. 
Clearly, there will be differences just as there are even in analytic resummation, where 
a much higher logarithmic accuracy can be achieved. Nevertheless, studying particular 
examples allows to quantify the deviations and convince oneself that they do not point 
to real errors either within the concept or the implementation. For this first study we 
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have chosen a process, which involves non-trivial colour exchange and massive partons, and 
requires cuts in order to even define the Born approximation cross section: top-quark pair 
production in association with at least one jet in hadron collisions at the Large Hadron 
Collider (LHC). There is one drawback of having top-quarks in the final state. A realistic 
description would require the inclusion of decays. We did not want to consider off-shell 
effects and the narrow width approximation was also out of the question for technical 
reasons. In such cases, one relies on the shower program to decay the quarks. This is a crude 
approximation, which neglects spin correlations. Unfortunately, DEDUCTOR does not even 
provide that. Thus, we are simulating a process with stable top-quarks in the final state. 
On the other hand, in our comparisons we have only taken the perturbative evolution into 
account. In other words, we have switched off the hadronization and multiple-interaction 
models available in common Monte Carlo generators (though not in DEDUCTOR). This 
allows us to quantify differences at the same level of approximation. 

We should point out that our matching implementation is restricted to leading colour. 
DEDUCTOR in the version we have used contains a more advanced colour approximation. 
Nevertheless, as long as full colour functionality is not available, we have decided to simplify 
our work. On the other hand, working at leading colour is appropriate for comparisons with 
other systems, which are also only certified to have this accuracy. 

Ours is the first phenomenological study performed using the Nagy-Soper shower, which 
does not involve the authors of the concept. We believe that this independence is important 
and actually proves that DEDUCTOR is ready to use by outsiders. It also gave us the 
motivation to prepare a summary of the main components of the Nagy-Soper parton shower 
including a list of ambiguities, which correspond to the places, where modifications might 
be expected once more experience in practical applications is gathered. We hope that this 
part of our paper will be useful to those, who do not necessarily want to read the hundreds 
of pages of the original publications to gain a basic understanding of the concept. 

The present text is organised as follows. We start by a summary of the Nagy-Soper 
parton shower concept. We are explicit as far as those parts are concerned, which have 
an impact on the construction of the matching procedure. We also point out some special 
features like the modified parton distribution function (PDF) evolution for example. The 
next Section deals with the matching itself. We consider both the case of processes, which 
have well-defined Born approximation total cross sections, and the case of processes, which 
require the specification of cuts already at this level. Subsequently, we describe the details 
of our implementation within HELAC-DlPOLES and the interface to DEDUCTOR. Finally, 
we show some results for ttj production at the LHC and compare them to results obtained 
with other systems. We close the main text with some short conclusions. An appendix 
contains a more thorough discussion of the matching in the leading colour approximation. 

2 A parton shower with quantum interference 

In this Section we review the Nagy-Soper parton shower concept, which allows for parton 
state evolution to include both spin and colour correlations. Neverthless, the main focus of 
our discussion will lie on the colour evolution. In fact, the exponentiation of non-diagonal 
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colour matrices is a long-standing problem. Notice that an independent attempt to include 
full colour evolution in a parton shower was already made in Ref. [11]. 

2.1 Quantum configuration space 

A generic 2 —>• m process is defined by two initial state partons a and b 1 and 
final state particles. Each particle is described by a set of quantum numbers to define the 
flavour fi, spin Sj and colour a of the particle and its momentum p t . Initial state partons 
are described by momentum fractions rj a and 775 . Thus, a complete m-parton ensemble, 
i.e. the probability distribution over possible quantum states, can be parametrized by the 
following set of quantum numbers 2 

{fh f 1 S) c}m = fa-i C a ], [%, /ft, Sft, Cft], [pi, /l, Si, Cl], ..., [p m , fmi Smi Cm]} . (2-1) 

Special emphasis should be given to initial state partons. In the Nagy-Soper shower, charm 
and bottom quarks are allowed to be massive. In the presence of parton masses, the initial 
state momenta are parametrized by the hadron momenta pa and pb, and the momentum 
fractions p a and % 
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m f 

Pa = VaPA H- —PB , (2.2) 

VaS 

2 

m f 

Pb = VbPB H- —PA , (2.3) 

VbS 

where p 2 A = p 2 B = 0 3 * and s = 2{pa • Pb)- In the following, we will describe the treatment 
of spin and colour in the parton shower. A matrix element A4 can be viewed as a vector 
f}m)) in colour <S> spin space and can be resolved into basis vectors |{s} m ) and 
|{c} m ) with complex expansion coefficients Ai({p, f, s, c} m ) 

I M({p,f}m))= M ({P’ IMm) ® |{c}m) 

{<s}m 

= Y - M (fr>/’' S > C }m) |{s,c} m ) . 

(s,c} m 

The basis in the spin space is orthonormal. Here, {s} m = {s a , Sb, s i,..., s m } represents the 
physical helicities of all particles. 

The colour basis is defined in terms of colour string configurations [36-40] , which allows 
for a straightforward connection with a hadronization model based on colour strings like 
the one described in Ref. [41]. In practice, a colour state |{c} m ) is represented by several 
strings {5i, ..., }, where Si can be an open or a closed string. An open string is denoted 

1 We treat the case of incoming hadrons. Colour neutral particles imply the usual simplifications, i.e. 
lack of parton distribution functions and no initial state evolution. 

2 The minus sign for the initial state flavour is only a convention, since all partons are considered as 
outgoing from the hard interaction. 

3 Even though parton masses are kept non-vanishing, one assumes that the proton mass is negligible with 

respect to the hadronic center-of-mass energy. 


(2.4) 
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by S = [Zi, I 2 , l n - 1 , ln \j where Zi is the label of a quark, l n of an anti-quark and I 2 , ..., l n -1 
are gluon labels. Furthermore, a closed string is given by S = (Zi,...,Z n ), where Zj labels 
only gluons this time. For brevity we will write h = i in the following. The label l of a 
parton refers to its colour index ap In case l labels a gluon, the colour index ai has values 
1,8. In case of quarks or anti-quarks a/ only takes the values 1,2,3. The set of colour 
indices { 01 , 02 ,-■■ , o n } is abbreviated in the following with {a} and allows to define the 
colour basis structures. 

The colour structure of an open string can be defined by 

fW(5 = [l,2,..,n—l,n]) = 1 [T a >T a * ■ • • T a ”~% ian , (2.5) 

Y N c Cp~ 2 

where T a is an SU (N c ) generator. Open strings are normalized to 

<S|S)^|tM (5) |2 = 1 . (26) 

{«} 

For closed strings we define 

= (1,..., n)) = —i=Tr[T ai r a2 , (2.7) 

\/Cp 

with normalization 

<S| S ),£|*M (5 )| 2 = 1 _(_±_)”-\ (2 .8) 

M 

A colour state |{c} m ) is represented by a product of these colour structures. Thus a state 
is normalized to 

({c} m |{c} m ) = n(5 fc |5 fc > . (2.9) 

k 

With these definitions the basis is not orthonormal in general. For instance 

({c'} m |{c} m ) = 5 (c} m ) + 0(l/N 2 c ) . (2.10) 

Therefore, the basis is only orthonormal in the leading colour approximation, N c —>• 00 . 

2.2 Quantum density matrix 

The basic object describing the parton shower evolution is the quantum density matrix p, 
which gives the “probability” 4 to find a certain parton ensemble {p, f} m . In this Section 
we give the relation of the quantum density to matrix elements. Using the notation of the 

4 Since p contains the full colour information, it can become negative for subleading colour configurations. 
Thus, strictly speaking, one cannot interpret p as a probability distribution. Nevertheless, the standard 
concepts from statistical mechanics apply here. 
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previous Section, one can write the expectation value for a completely inclusive observable 
F , as 


"1*1 - E ^ / Wl'>»> wi* /}-)> 

= ^2~i f [d{pJ}m]Tr[p{{P, f}m)F({P, f}m)} , 

'' ml 

m J 

( 2 . 11 ) 


where the sum runs over all final state multiplicities. Here, [d{p, f} m ] is the sum of all 
m-particle phase space measures for different flavour sequences {/} m . The factor 1/m! is 
necessary to account for identical contributions. The parton density functions evaluated 
at the momentum fraction rj and factorization scale p 2 F are denoted by f a /b(jh P F )- We 
average over initial state spins (factor 4 in the denominator) and colour (factor n c (i), with 
n c {q ) = 3 for quarks and n c {g ) = 8 for gluons). The trace in the second line corresponds to 
a sum over indices in the colour <S> spin space. The quantum density p is, therefore, given 

by 


p({pj}m) 




fa{Pa,P 2 F )fb(.Vb,PF) 

4 n c (a)n c (b) x flux 


EE \{S, C Jm ) P({P, /, S', C, S, C} m ) ({s', C'U , 

S,c s',c' 


( 2 . 12 ) 


where we used the expansion of | Ai({p, /} m )) according to Eq. (2.4). The quantum density 
matrix is a projector onto the different helicity and gauge invariant subamplitudes of the 
full quantum amplitude 


P({p,f, s',c', 


s,c} m ) = M*({p,f,s', 


c'}m)M({p, f,S,c}m ) 


fa(ria,P F )fb(.r]b,P F ) 
4 n c (a)n c (b) x flux 


(2.13) 


Using the quantum density matrix, the expectation value of the observable F is given by 


5 


V[F 1 = f’ S ’’ C ‘ S) f}m) ({«') c} m ) p({p, /, s' , c', S, c} m ) , 


(2.14) 

where the integration measure has been extended by the inclusion of the summation over 
spin and colour. It is useful to define basis vectors or statistical states \{p, f, s', c', s, c} m ) 
(see Ref. [16]), such that 


P({P, /) s', c', s, c} m ) = ({p, f, s', c', s, c} m \p) . (2.15) 

Defining additionally an abstract state (T| according to 

(F\{P, /) s', c', s, c} m ) = F({p, f} m ) ({s', c'} m |{s, c} m ) (2.16) 

5 We assume that F is a unit operator in colour <X> spin space. 
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and using the decomposition of the identity 


1 = 


Y2 “T / i d ({P’f’ S 'i C 'i S ’ C }rn)}\{pJ,s', c ',S,c} m )({p,f,s',c',S,c} m \ , (2.17) 


one can reduce Eq. (2.14) to a scalar product 

a[F\ = ( F\p ) . 

Finally, we define the total cross section measurement function (1| as 
{ l \{ p , f , s ', d , S , c } m ) = ({s'} m |{s} m ) {{ d } m \{ c } m ) 


(2.18) 

(2.19) 


2.3 Evolution equation 

The evolution equation describes the propagation of the quantum density matrix p from 
some initial shower time to to some final time tp. The initial time corresponds to the hard 
interaction, while the final to the formation of hadrons. Therefore, tp characterizes the 
physical scale at which parton emissions cannot be described perturbatively. There is no 
unique definition of shower time t as explained in more detail in the upcoming Section 2.7. 

The perturbative evolution is described by an operator U(tF,to). The expectation 
value of the observable F, including shower effects, is 

a[F] = ( F\p(t F )) = (F\U(t F ,t 0 )\p(t 0 )) . (2.20) 

U ( tp , to) is assumed to be unitary 6 in the sense that the total cross section ap is not affected 
by evolution 

°T = (1| p(t F )) = (l\U(tF,t Q )\p(t 0 )) = (l\p{t 0 )) ■ ( 2 -21) 

The evolution operator is the solution of the equation 

= mt) - V(t)]U(t, t 0 ) . (2.22) 

Here, 7ij(t) describes the transformation of a state {p, /, s', c', s, c} m to another state 
{p, /, s', c', s, c} m +i by the emission of a resolved particle. This transition between the 
states is of course constrained by overall momentum conservation. On the other hand, 
V(t) describes the unresolved/virtual emission and, therefore, does not alter momentum or 
flavour configurations of the particles. Nevertheless, it can change colour configurations, 
which can have an effect on further emissions. It can be further decomposed into a colour 
diagonal, Vp(t), and a colour off-diagonal part, Vs(t) 

V(t) = V E (t) + V s (t) . (2.23) 


Traditional parton showers correspond to the large N c limit. The colour structure is, 
therefore, always diagonal (Vs(t) —>■ 0). In that case, Eq. (2.22) can be solved as 

U(t, t 0 ) = N(t , to) + f dr U(t, t)'H i (t)N(t, t 0 ) , (2.24) 

Jto 

6 A non-unitary evolution can be used to resum soft gluon effects. 
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where N(t,to) is the Sudakov form factor (a number) defined as 


N(t, t 0 ) = exp 



(2.25) 


In case of a non-trivial colour evolution, the exponentiation of a non-diagonal matrix 
is cumbersome. Instead, only the colour diagonal part, V^;(t), is exponentiated while the 
off-diagonal part, Vs{t) is treated perturbatively on the same footing as Ri(t) 


U(t, t 0 ) = N(t , t 0 ) + f dr U(t , r) [W/(r) - V s (r)] N(t, t 0 ) , (2.26) 

J to 

with ^ 

N(t, to) = exp ^ dr V E (r)^ . (2.27) 

This is the final evolution equation of the Nagy-Soper parton shower. Notice, however, 
that the current implementation in DEDUCTOR involves some approximations which are 
documented in Section 2.6. 

In the following we provide the definitions of the real splitting operator "H/(t) and 
the virtual operator V(t). As we will see, the unitarity condition on U(t,tp) will allow to 
determine V(t) in terms of T-Lj(t). 


2.4 T~Li(t) - Real splitting operator 

The real splitting operator describes the transition of an m-particle ensemble to an (m+1)- 
particle one, by splitting a particle into two, 

{p, f, s ', c', s, cjm {p, /, s', c’, s, c} m+i . (2.28) 


This splitting is constrained by flavour and momentum conservation. The parameters of 
the transition require further specifications, because they depend on a particular choice of 
momentum mappings, shower time and splitting functions. The behaviour of the splitting 
operator is fixed in the infrared regime in order to reproduce the correct singular limits of 
QCD amplitudes. 

After emitting a particle, one has to modify the momenta in the event to preserve 
momentum conservation and the on-shellness of all particles. This is achieved by momentum 
mappings Ri, which require three additional variables parametrizing the momentum of the 
emitted particle, T/, and a variable, \h specifying its flavour f m+ 1 . Thus, the new momenta 
and flavours are given by 

{P, f}m+ 1 = Rl({p, f}m , r l,Xl) , ( 2 - 2 9) 

where l € {o, b, 1 , ...,m}. The inverse of this transformation is denoted by Qi 

Qi{{p, f } m+i) = {P, f}m ■ (2.30) 

The Nagy-Soper parton shower makes use of a global momentum mapping [16], i.e. all 
final state particles’ momenta are modified to account for the momentum of the emitted 


particle. The freedom inherent in the definition of Ri may be exploited to improve the 
resummation for certain observables. For instance, the study presented in Ref. [19] has 
shown, that the pt spectrum for Drell-Yan Z-boson production depends strongly on the 
momentum mapping for initial state parton splittings. The original Ri from Ref. [16] has 
been found inadequate and subsequently modified in Ref. [22], The modified version has 
been adopted in Deductor. 

One can define a momentum mapping operator Vi satisfying two conditions, the first 
one being 

—7 / [d{p',f'}m\ ({pJ}m+l\'Pl\{p',f'}m)g({p',f'}m) = g({p, f}m) , (2.31) 

m\ J 

where g({p, f}\) is an arbitrary test function and {p, f} m is determined by Eq. (2.30). The 
second condition is 

7 - 3 -jTT [[d{pj}m+l\g({pj}m+l) ({p, f}m+l\Vi\{p, f}m) 

[m + 1 )! J 

1 ^ f , (2.32) 

= (^TI)_E JrtuiMf'}«+!), 

fm+lC-Xl 

where {p' , f'}m- 1-1 is determined by Eq. (2.29). Vi contains a 5-function, which links the 
(m + l)-particle kinematics to the m-particle kinematics with R\. The two conditions imply 
that the Jacobians resulting from the integration over the 5-function in Eqs. (2.31) and 
(2.32) have been absorbed in Vi- 

The evolution of the density matrix can be determined by studying the factorization 
of QCD amplitudes in the soft and collinear limits. We remind that in the limit when two 
partons become collinear, pi || p m+ i, the amplitude factorizes as 

\M({p,fU +1 )) « T }(/, h + f m+ i)V?({pJ} m+ i) | M({p,f}m)) , (2.33) 

where ?)'(_/) —> fi + f m +i) is an operator in colour space (see Ref. [16]). V t '({p, /} TO +i) 
is the splitting operator in spin space. Contrary to traditional parton showers, the Nagy- 
Soper formulation does not use the Altarelli-Parisi splitting kernels [42]. Instead, the re¬ 
quired functions are derived directly from matrix elements without taking the collinear 
limit first [16]. This allows, for example, for a direct access to polarization information. 
Furthermore, the functions are not singular in the soft limit [17]. In the soft limit, on the 
other hand, when the momentum of a gluon vanishes, a similar approximation is valid 

\M({p, f} m+ 1 )) « ^ t/(/, ^h + f m+ i)vr h \{p, f}m+ 1 ) I M({ p , f } m )) , (2.34) 

/ 

where l runs over all partons and V) soft is an eikonal factor. 

Combining these approximations, we may write the following evolution equation for 
the quantum density matrix 

|Pm+i) = I Pm) ■ (2.35) 

l 
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Consequently 7~Li(t) is defined by the splitting operators Si at a fixed shower time Ti({p, f } m ) 
(see Section 2.7) 

Ui(t) = ^> - 7K(P> /}m)) • (2.36) 

i 

We shall not reproduce the exact form of 5/ in the most general case, as it can be 
found in Ref. [16]. Nevertheless, we would like to point out that there is an ambiguity in 
distributing the soft limit among the contributions with different momentum mappings. In 
order to illustrate this issue, we give the form of 'Hz(t) in the slightly less involved spin 
averaged case 


{{P, f, c , c} m+ i\Hi(t)\{p, /, c', c} m ) 

= {m + l) ^ 2 6 (t-Tl({p,f}m )) ({pJ}m+l\Vl\{p, f}r 


l,k 


n c (a)n c {b)rj a rj b fa/A^Pa, Mf) ^f) 

n c (a)n c (b)f) a fi b f a /A(Pa, PF)fb/B(.Vb, Pf) 

l r. 


X - 


2 ^ki(l-S frn+ iJwii({p,f} m+1 ) 


+ SklS 


/ ra +1 


Wll({p, f}m+ 1) - wft konal ({p, /} m+1 ) 


— (1 - Ski)S^ 
X 


ra+ 1 ,9 


Alk({p}m+l)wf k POle ({p, f} m+ 1 ) 


({c 7 , c} m _|_i|T^(/( -> fi + /m+l) <8> T k {fk fk + /m+l)|{c / , c} m ) 
+ ({c / , c} m +l|rj(/fc — > fk + /m+l) SS T[(fi — > fl + / m +l)|{c / , c} m ) 


(2.37) 


where the sum runs over two partons l and k, not necessarily different, and l is the emitter. 
The function wu is the spin averaged squared splitting function, while yj^J- konal i s its soft 
approximation (see Refs. [17, 21]). Furthermore, 


wt P ° le ({pJ}m+ 1 ) = -4vra s 


{{pm+ 1 • Pl)Pk ~ {Pm+ 1 • Pk)Pl) 2 

{Pm +1 •p0 2 (Pm+l • Pk ) 2 


(2.38) 


which, for massless l and k , reduces to the well known squared eikonal factor [43]. The func¬ 
tion Aik{{p] m +i) removes the singularity when the partons m + l and k become collinear. 


A-lk({p}m+ 1 ) 


_ (Pm+l • Pk){Pl • Q) _ 

(Pm +1 • Pk)(Pl • Q ) + (Pm+l • Pz)(Pfc • 0) ’ 


(2.39) 


where Q is the total final state momentum. The ambiguity we have mentioned is due to 
the fact that the product Aikwfl P ° le must only satisfy two conditions: 1) it may not have a 
singularity in the collinear limit p m +i || Pk : 2) it should correctly reproduce the singularity 
at Pm+l —> 0. Formulae (2.38) and (2.39) only correspond to one specific choice. 
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2.5 V(t) - Virtual splitting operator 

As advertised, the virtual splitting operator, V(t), can be determined to a large extent from 
the real splitting operator fii{t) by the requirement of unitarity 

(l\n I (t)-V(t)\p(t))=0, (2.40) 

which should be valid for any |/?(t)). This can be interpreted with the help of the optical 
theorem. For example, we can schematically write 



-V-' 

(l|ft/(t)|p) (l|V(t)|p) 


Notice that V(t) does not contain the exact virtual corrections. Furthermore, it turns out 
that /, s', c', s, c} m ) does not contain spin correlations due to the integration 

over the azimuthal angle of the emitted partons [16]. This allows us to write 

(l\Hi{t)\{p, /, s', c', s, c} m ) = 2 ({s'} m |{s}m) ({c'} m \h(t, {p, f} m )\{c}m) . (2.42) 


In consequence, the freedom in defining V(t) is restricted to the colour structure. In princi¬ 
ple, a virtual correction may be applied either on the c' or c indices of the | {d, c } m ) state, 
see Fig. 1. This leads to the following 





Figure 1. Modification on the colour structure by the inclusion of the virtual operator V(t). 


V(t) = (h + i(j>) <S> 1 + 1 (hf — i(j>) , (2.43) 

where 0 is a colour dependent phase, which can not be determined from real radiation 
corrections. Its presence is a consequence of a Coulomb gluon exchange [21] 

I 2 2 

<Kt, {p, f}m) = -27T ^ ^ — T k ■ T; , v kl = J 1 - , (2.44) 

^47 TVkl V [Pk-Pl) 

where the sum runs over l and k, both either in the initial or final state. Examples of the 
relationship between 'Hj(t) and V(t) are depicted in Fig. 2. The exact specification of V(t) 
can be found in Refs. [16, 21]. 
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Figure 2. Correspondence between a real radiation contribution as given by Hi{t) and its counter¬ 
part in V(t), for a direct splitting (upper line) and an interference contribution (lower line). In both 
cases the left part of the first graph corresponds to l = 1. For the right part k = 1 (upper graph) 
or k = 2 (lower graph). In the case of the interference contribution, by convention, we attach the 
virtual correction to the amplitude of the parton k. 


2.6 Colour and spin evolution 

An essential part of the action of is the change of the colour state. Starting from the 

two types of colour strings Eqs. (2.5) and (2.7), the application of T} ifi fi + f m + 1 ) on 
the c' index and of Tk(fk —>■ fk + fm+i) on the c index of the |{ c',c} m ) state results in a 
linear combination of the same basis vectors. This follows from the relations 


ifabcT c = 


T J CL rpd _ 

ij 1 kl — 


j _ 'j^drj^b _ rj-ibrj-ia 


k jy &ij &kl 


(2.45) 

(2.46) 


which can be translated into similar relations for T^\fi —>■/; + / m +i), see Ref. [16]. 

In principle we have to consider diagonal, {d} m = {c} m , and off-diagonal, { c'} m 7 ^ 
{c} m , colour configurations, see Fig. 3 for illustration. Notice that the colour diagonal con¬ 
tributions contain the leading colour contribution, while colour off-diagonal contributions 
are subleading. The latter can be identified by crossing of gluon lines. Traditional parton 
showers are restricted to the leading colour approximation (or at best to colour diagonal 
contributions). On the other hand, the Nagy-Soper parton shower allows for the evolution 
of subleading colour configurations. In Figs. 4, 5, 6 , 7, and 8 , we give several examples 
of colour evolution induced by direct splitting (application of Tj <£> T{) and interference 
contributions (application of T) ® T)., where l 7^ k). 

The solution of the evolution equation, Eq. (2.26), requires a decomposition of the 
virtual operator into a diagonal and an off-diagonal colour part. This can be achieved 
with the help of the LC+ approximation introduced in Ref. [21]. The main features of this 
approximation are 


1. Exact colour treatment of the collinear and soft-collinear limits; 
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Figure 3. Example colour configurations for the q(a)q(b)g(l)g(2) state. 
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Figure 4. Example of a direct splitting g —> gg starting from a diagonal colour configuration. The 
left part of the first graph corresponds to l = 1, while the right part to k = 1. Moreover, m + 1 = 3. 



Figure 5. Example of a direct splitting g —> qq starting from a diagonal colour configuration. The 
left part of the first graph corresponds to l = 1, while the right part to k = 1. Moreover, m+ 1 = 3. 


2. Leading colour approximation for the pure soft limit (a part of the subleading contri¬ 
butions is kept nevertheless). 

According to the previous Section, the virtual splitting operator is defined in terms 
of the real splitting operator. For the latter, the LC+ approximation amounts to the 
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Figure 6. Example of a direct splitting g -+ gg starting from an off-diagonal colour configuration. 
The left part of the first graph corresponds to l = 2, while the right part to k = 2. Moreover, 
m + 1 = 3. 



a a 

1 1 

mm' 2 3 * otToo 
cro'SC' ,J 2 'cro'cro 

v_ 6 b 





'O'O'O'O 

'"O'O'O'O" 


LC+ 

—» 



7 HRHT 2 2 

V_ 66 _^ 



"oootr 2 3 TJ5J5C- 

TKRHT 2 2 ^CRJOTJ- 

v_6 b _y 


Figure 7. Example of an interference contribution of two g -+ gg splittings starting from a diagonal 
colour configuration. The last line represents the LC+ approximation of the original colour struc¬ 
ture. The left part of the first graph corresponds to l = 2, while the right part to k = 1. Moreover, 
m + 1 = 3. 


introduction of a projector C(l,m+ 1) acting on the colour states |{c} m +i) in Eq. (2.37) as 

— * fk + fm+l) |{c}m) —> c(l, m + l)7fc(/fc —t fk + fm+l) |{c}m) , (2.47) 

({c^ml Tk{fk fk + fm+l ) ~^ ({ c/ }m| Tk{fk fk + /m+l)C^(|, m + 1) , (2.48) 
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Figure 8 . Example of interference contributions of g —» gg and q —^ qg splittings starting from 
a diagonal colour configuration. Upper graph: the left part of the first graph corresponds to l = 2, 
while the right part to k = b. Lower graph: the left part of the first graph corresponds to l = a, 
while the right part to k = 1 . Moreover, m + 1 = 3 in both cases. The last line represents the LC+ 
approximation of the original colour structure of the lower graph. 


where 


C(l,m+ 1) |{c} m+ i) 


|{c}m+l) 
11 |{ c }m+l) 
0 


partons l and m + 1 form a qq pair , 
partons l and m + 1 are colour connected , 
otherwise . 


(2.49) 


Colour connected partons are partons, which are direct neighbours on a colour string as de¬ 
fined in Section 2.1. The definition Eq. (2.49) guarantees that direct splittings are treated 
exactly. This is also true for the evolution of off-diagonal colour configurations, as de¬ 
picted in Fig. 6. In other cases, a part of the contributions is removed, as can be seen in 
Figs. 7 and 8. 

The colour projector C(l,m + 1) defines a decomposition 

Ui{t) = 'H 1 f c+ {t) + AH^t) . (2.50) 

All corrections to the LC+ approximation are of order 0(1/N (?) and are treated perturba- 
tively. It turns out that the virtual splitting operator V LC+ (t ) defined in terms of TL I f C+ (t) 
is colour diagonal [21]. Therefore, we define 

V E (t ) = V LC+ (t) , V s (t) = AV(t) = V(t) - V LC+ (t) . (2.51) 


We can also define an approximate evolution operator 

U LC+ (t, t 0 ) = N LC+ (t, t 0 ) +[ dr U LC+ (t, t)TL I f c+ (t)N lc+ (t, t 0 ) , (2.52) 

Jt 0 
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where ^ 

N LC+ (t,to) = exp ^ dr V LC+ (r)^ . (2.53) 

On the other hand, the exact evolution operator is 

U(t, t 0 ) = U LC+ (t, to) + r dr U{t , r) [Aftj(r) - AV(r)] U LC+ (r, t 0 ) . (2.54) 

Jt 0 

The practical solution of these equations is described in Ref. [21]. The current implemen¬ 
tation of Deductor relies on Eq. (2.52). 

Let us comment on the logarithmic accuracy of the Nagy-Soper parton shower. Con¬ 
sider an observable O , given by the expansion 

(O) = Y c(n, 2 n)a n s L 2n + Y c(n, 2n - 1 )a n s L 2n ~ 1 + • • ■ , (2.55) 

n n 

where L denotes a potentially large logarithm, e.g. L = log (s/p^), where s is the scale 
of the hard process, and px is the transverse momentum of some final state, for example 
a pair of oppositely charged leptons from the Drell-Yan process. The evolution operator 
U(t,to) reproduces the coefficients c(n, 2 n) and c(n, 2 n— 1 ) without any approximations in 
colour. On the other hand, U LC + (f, to) misses part of the c(n,2n — 1) coefficient due to 
subleading colour contributions from wide angle soft radiation. For this reason, the LC+ 
approximation is only exact at leading logarithm (LL), see Fig. 9. Inserting a factor of 
[A 7~Li(t) — AV(t)] generates the remaining contribution to the coefficient c(n, 2n — 1). In 
consequence, the evolution operator 

U(t, to) = U LC+ (t , to) + f dr U LC+ (t, t) [Aftj(r) - AV(r)] U LC+ {r , t 0 ) , (2.56) 

■Jto 

is accurate at next-to-leading logarithm (NLL). 

As pointed out before (see discussion around Eq. (2.42)), the virtual splitting operator 
is diagonal in spin. Therefore, spin dependent effects are only induced by the action of 
7ij(t). It has been proposed in Ref. [18] to include them on top of the spin averaged 
shower. To this end, we note that the necessary modification in 7ij(t), Eq. (2.37), amounts 
to the replacement 

®lk({p,f}m+ 1 ) -» ({s',s}m+l\Yl k ({p,f}m+l)\{s',s} m ) $lk({P,f}m+l) , (2.57) 


where 


*&lk({p, /}m+l) 


and Y lk ({pJ} m + 1 ) is just 
functions, see Ref. [18]. 


= 5 ki (1 - 8j m+i<g ) wu({pj } m+ 1 


+ $kl 8 


fm+119 


Wll{{P,f}m +1 - W 


■eikonal 

ll 


(2.58) 


- (1 - 6 k i) <5; m+li9 A lk ({p m+1 }) wf l k pole ({pJ} m+ 1 ) , 

the ratio between spin dependent and spin averaged splitting 
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LL NLL 

Figure 9. Illustration of the logarithm counting. One step in the vertical direction is given by an 
insertion ofU LC+ and a diagonal step is given by an insertion of AHi(t) — AV(t). One can see, 
that two insertions of AHi(t) — A V(t) only contributes to the coefficient of a" L 2n ~ 2 . 


The parton shower evolution starts with a quantum density matrix 

(-fa, /, s', d, s, c} m \p(t 0 )) = (fa, /, c', c} m \p(t 0 )) (g) ({s', s} m \p s ({p, f, d, c} m )) , (2.59) 

with a factorized spin dependence, |p s ), with respect to the spin averaged quantum density 
\p(to)). The latter subsequently undergoes spin averaged evolution. On the other hand, 
after (n — m) real emission steps 

I P TL ) 1 ( 1 f { U,') 2 (fa, f } ,km+i (fa, /} m+1 )|P (fa, /, c , c }m)) • 

(2.60) 

Assuming that the hnal observable does not depend on the spin configuration, we obtain 
°[F] = X! TT / l d iP’ f’ c 'i c }*\( F \{P’ f’d,c}\)({p, f,d,c}x\p(tF))(hpin\P\) , (2.61) 

A ' J 

where 

(hpin\p{)= ^2 (lspin\{s',s}x)({s',s}x\p{) = ■ ( 2 ’ 62 ) 

(s',s} A {s} A 

2.7 Shower time 

The forward evolution in the shower time is used to enforce an ordering in some chosen 
kinematic variable. This is necessary in order to correctly resum leading logarithms of 
infrared sensitive quantities. In traditional parton shower programs the following ordering 
variables are used: 

virtuality: the invariant mass of two daughter partons produced through the splitting (e.g. 
Pythia6Q [44]), 
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transverse momentum: the transverse momentum of the daughter partons with respect 
to the mother parton (e.g. PYTHIA8 [45]), 

angle: the angle between the momenta of the daughter partons (e.g. IiERWIG [46, 47]). 


The construction of a new parton shower allows to re-evaluate the advantages of these 
variables, and replace them if necessary. In Ref. [22], the following variable has been 
proposed 


A i = 


±p m+ 1) 2 - rrvf 


■Qo 


(2.63) 


2 Pi ■ Qo 

where pi is the emitter momentum after emission (daughter parton), p m+ i the emitted 
parton momentum (daughter parton), pi the emitter momentum before emission (mother 
parton with mass mi) and Q q is the total final state momentum. The minus sign between 
pi and Pm+i in Eq. (2.63) applies to an initial state splitting. The shower time with Af 
ordering is then given by 


7K(p> f}m) = - log ■ (2.64) 

The particular form of Af may be justified as follows [22]. Consider a splitting of a 
final state parton with momentum po into two partons with momenta p\ and P 2 ■ After 
splitting the mother parton has non-vanishing additional virtuality Vq (= Pq — mg), while 
both daughter partons are on-shell. In subsequent steps, both daughter partons will split 
as well, and acquire additional virtualities vf and (= pf — rnf). Clearly, Vq will be 
modified by this procedure. The validity of the on-shell approximation used in the first 
step necessitates the effect on Ug to be negligible. This can be translated into the following 
conditions [22] 


vl 


< 


2pi • Qo 2p 0 • Qo 


and 


V 2 

v 2 


< 


2 P 2 ■ Qo 2po • Qo 


(2.65) 


Up to normalization, the terms 2p"-Q 0 are ec l ua l to A^. A similar discussion holds for the 
initial state shower. We note that one of the consequences of Af ordering is an enlarged 
phase space for initial state splittings compared to px ordering. 


2.8 Consequences of backward evolution 

Initial state splittings are generated by the backward evolution formalism first introduced 
by Sjostrand in Refs. [2, 3]. In this approach, the Sudakov form factor contains a ratio 
of parton distribution functions to guide the evolution of the momentum fraction during 
splitting. Nevertheless, as long as mass effects are neglected and the splitting functions used 
in the shower and PDF evolution are the same, one can factor out the non-perturbative part 
completely. This matches the treatment of final state radiation, where hadronization only 
happens after showering. Unfortunately, the splitting functions used in the Nagy-Soper 
parton shower are very different from those of Altarelli-Parisi away from the collinear limit. 
Thus, non-perturbative information will inevitably affect the evolution. 

In the presence of mass effects, however, the non-perturbative input can not be factor¬ 
ized, independently of the approach to the parton shower (traditional or otherwise), unless 
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the evolution of PDFs is modified. We describe the required modification after Ref. [23], 
which has been implemented in Deductor. 

In principle we would like to decompose the density matrix as 


Ip) — £(^)|Ppert) i 


( 2 . 66 ) 


where |p pe rt) would contain no non-perturbative information if splittings were restricted to 
the quasi-collinear limit. Here, the luminosity operator C(t) is 


C{t )\$ m ) 


fa/A{Va^A e *) fb/B^b, P? B e *) 

An c (a)n c (b)4r] a r] b p A ■ p B 


(2.67) 


where we used a shorthand notation for the statistical states, = {p, f, s', c', s,c}\. The 
scale p‘f is given by fij = 2 pi ■ Q o, where p % are the hadron momenta and Q o the total hnal 
state momentum. Therefore, the PDFs are evaluated at the scale /r 2 e ~ 4 = \(pi~p m +i ) 2 ~Jnf | 
which is the virtuality of the splitting. Let us define a perturbative real splitting operator, 

^pcrt(^), t ]- lr{)U g} 1 

Mi(t) = C(t)U¥ irt (t)£- 1 (t) . (2.68) 

Unfortunately, V(t) contains a convolution of splitting functions with PDFs, because it is 
an integral of Eq. (2.68). It is thus impossible to trivially factor out the non-perturbative 
contribution. Nevertheless, we can define a quasi-perturbative virtual splitting operator, 
V pert (i), via the evolution equation. Indeed, substituting Eq. (2.66) in Eq. (2.22) we obtain 


IPpert) + ■£(£) ^ ^IPpert)^ — [Hl(t) V(t)] C(t) |p P ert) ■ (2.69) 

This allows us to write 


= [«?“(*) - v-“(t)] | 


(2.70) 


with 

V pert (f) = V(t) + £” 1 w(|^)) , (2-71) 

where we used the fact that C(t) commutes with V(t) since it does not change the momenta. 

To ensure that V pert (t) is independent of PDFs in the quasi-collinear limit, the splitting 
functions obtained by deriving C(t) must match the quasi-collinear limit of V(t). The 
respective functions are given in Ref. [23]. As expected, they only differ from the Altarelli- 
Parisi kernels at non-zero mass. They also imply a modified evolution of the PDFs. The 
effect is mostly visible on the 6 -quark PDF. Clearly, only processes predominantly generated 
from initial state 6 -quarks will be influenced. 

Modifying the evolution of parton distribution functions may be used for other purposes 

( 2 ) 

as well. For instance, adding a higher order term P- as follows 


dfa/A{Va,^ 2 ) 

dlog(/r 2 ) 


Y J ^7 aS ^J^ P a a(z, P 2 /Z) fa /A {ri a /Z , p 2 ) 


+ 


^2 J — 


P !S( z ^ 2 / z )fa/A{Va/z,P- 2 ) , 


(2.72) 
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with 


p{ ah^ Z ^ 2 / Z ) = - 27 r /5olog(A R)Paa{z,H 2 /z) , 

where /3q = (33 — 2ny)/(127r) and 


A r = exp - 


Ca {67 — 37t 2 ) — 10 ra/A 

3(33-2 n f ) ) 


(2.73) 


(2.74) 


where rif is the number of light active flavours, allows to incorporate the cusp anomalous 
dimension at NLL for processes with less then three coloured particles [48]. 


2.9 A summary of ambiguities 

Let us conclude this Section with a list of ambiguous components of the Nagy-Soper parton 
shower: 


1. Momentum Mappings 

2. Splitting functions 

3. Soft partition function 

4. Colour treatment 

5. Spin treatment 

6 . Shower time 

7. PDF evolution 

We stress that this shower is a rather recent construction still under development. It is to 
be expected that there will be changes to any of the particular solutions discussed in the 
previous Subsections. 


3 Matching NLO matrix elements to the parton shower 

Matching NLO calculations with parton showers is a widely explored subject and there exist 
several matching schemes, the most popular being POWHEG [28, 29] and Mc@NLO [26, 27] 
(for other proposals see e.g. [13, 49-56]). A general comparison between these two can be 
found in e.g. [57]. In order to benefit from the recently implemented subtraction scheme 
based on the Nagy-Soper parton shower splitting kernels [33], we chose the Mc@NLO 
formalism. Notice, that the implementation of Mc@NLO in SHERPA [57] also benefits from 
the consistency between the parton shower (Catani-Seymour) and subtraction terms in the 
matching. However, once the implementation of the Nagy-Soper parton shower will be 
made accurate at subleading terms in the colour expansion, our matching procedure will 
allow for a similar accuracy in the matched sample. This may also be possible in the future 
within SHERPA, see Ref. [11], Before we discuss the challenges of parton shower matching, 
we give a brief overview of the general objectives. 
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Preserving the NLO cross section normalization: When considering an inclusive 
observable F, the fixed order normalization of the cross section should be preserved. 
We thus require 

(F\U(t F ,t 0 )\p(t 0 )) = a NLO [F] . (3.1) 

Defining event samples at NLO level: Matching to parton shower is the only way 
to define events at NLO. As we will see, matching a la Mc@Nlo imposes a shower 
specific subtraction scheme. Without matching, the weights of the real matrix element 
and the subtraction terms are not coupled kinematically unless in a strict limit, and 
diverge separately. Due to the matching scheme those weights are combined and one 
obtains real emission phase space configurations with a finite, but not necessarily 
positive, weight. The matching renders the virtual corrections finite as well. 

Next-to-Leading Logarithmic accuracy for infrared sensitive observables: In¬ 
frared sensitive observables, which are affected by large logarithms, L, of some kine¬ 
matic variable at fixed-order, are replaced by resummed predictions. Traditional 
parton showers only resum the leading logarithms (LL) a™L 2n , whereas the Nagy- 
Soper shower also allows to resum the next-to-leading logarithms (NLL) cd?L 2n_1 . 
Therefore, we require differential distributions to be accurate at the NLO+NLL level. 

High pt emissions according to matrix elements: The parton shower is valid in 
the soft and collinear regimes, the description of high px emissions is, therefore, not 
reliable. On the other hand, NLO calculations are valid in this region. It is desirable 
to recover the NLO prediction for high px emissions despite applying the parton 
shower. 

Connection to low energy physics: The parton evolution down to a low scale allows to 
include non-perturbative phenomena due to hadronization and multiple interactions. 


3.1 The quantum density matrix at next-to-leading order 

For a generic 2 —>• m process at NLO, one can write the quantum density matrix in a 
perturbative expansion in a Sl according to 


Ip) 


\^) 

Born, 0(1) 



Virtual, 0(a 3 ) Real, 0(a 3 ) 


+@( a s) ■ 


(3.2) 


Note that we count the leading order contribution as order 1 in the strong coupling a s . 
| pm 1 ) and j) are tree level matrix elements, whereas \pm) is the one-loop amplitude. 
The definitions of these densities are analogous to the definition given in Eq. (2.12). So far 
Eq. (3.2) suffers from infrared divergences, however, the expansion of \p) in a s is defined 
within dimensional regularization. Based on this quantum density matrix, the expectation 
value of the observable F including shower effects reads 


^ i r 

a[F] PS = (F\U(t F ,t 0 )\p) = E a? / [d*xm$ x ){§ x \U(t F ,t 0 )\p) , 

X=m 


(3.3) 


where = {p, f,s',c',s,c}\. The quantum density | p) accounts for the hard matrix 
elements for A = m and A = m + 1. This naive description of the cross section suffers 
from double counting, as we will show using the iterative solution to the evolution equation 
expanded to 0(a s ). Indeed, there is 

\p(t F )) = U(t F ,t 0 )\p)fa\p)+[ dr — V(t)] \p$) + O(a^) . (3.4) 

Jt 0 

As we can see from the unitarity condition (1| — V(r)] = 0, the total cross section 

(l|/o(ijr)) is conserved. On the other hand, one does not recover the NLO prediction for 
inclusive observables F because in general (F\ — V(r)] ^ 0. Even without this 

requirement the result is not correct, since it contains the first emission contributions twice, 
once from the real emission quantum density \pf n > + ]), and once from the parton shower 
approximation 'Hi(r)\p'm)- We will now show how to overcome this problem. 


3.2 Matching fully inclusive processes to parton shower 

Let us begin with a class of processes, which have a well defined total cross section at leading 
order, e.g. pp —>■ tt or pp —>• W + W~. Working along the line of the Mc@Nlo scheme, 
we notice that the additional parton shower contribution in Eq. (3.4) can be cancelled by 
including appropriate counterterms. Let us, therefore, introduce the following modified 
quantum density matrix 

Ip) = \p)~ f dr [^/( T ) - V ( r )l I Pm) + °( a2 s ) • ( 3 - 5 ) 

Jt 0 

The total cross section (l|p) = (l|p) = a NLO is unchanged, due to the unitarity condition. 
On the other hand, considering U(tF,to)\p) and expanding the evolution equation shows 
that the undesired parton shower contributions are cancelled up to 0(a s ). This cancellation 
is non-trivial, because the modified quantum density matrix |p), now depends explicitly on 
the parton shower splitting kernels, the momentum mappings, the ordering of the emissions 
and the choice of the starting shower time to- 
For an infrared safe observable F, we have 


v[F] = —. f [d$ m }(F\U(t F ,to)\$m) \($m\ P W) + ($ m |^ } ) + [ tF dr($ m |V(r)|pW) 
m - J L J to 


+ 


1 


(m + 1)! 


[d3>m+l] {F\ U(tF, to) 14>, n+ i) 


(^m+i|pJ2hi) - [ dr($ m+ i|^/(r)|p^ ) ) 

Jt 0 

(3.6) 


In the Mc@NLO approach, the parton shower splitting kernels are used to provide sub¬ 
traction terms for the infrared divergences in |pin ) and ). Thus, the infrared cutoff 

tp may be removed by taking the limit tp —>• oo. Even though this is a source of a mis¬ 
match between the fixed order and the shower calculation, it is numerically small due to 
the exponential damping by the Sudakov form factor as discussed in Ref. [2]. In the real 
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subtracted cross section, described in the second line of Eq. (3.6), we use the definition of 
the real splitting operator from Eq. (2.36) and write 


POO _ POO 

/ dr Ui{t) = y~]S; / dr 6 (t - f/)0(r - t 0 ) = V] S/0(t; - t 0 ) . 
Jt 0 , Jo , 


(3.7) 


Here the sum runs over all external legs and S i is the total splitting kernel for a given 
external leg l. We want to emphasize that S/ also contains non-singular contributions like 
the massive g —?• QQ splitting. Therefore, the operator T~Li{t) contains more than is needed 
for a subtraction scheme. The parameter f/ is the shower time defined in Section 2.7. 
Hence, 0(f; — to) represents the ordering of the emissions and the to dependence introduces 
a dynamical restriction of the subtraction phase space. The choice of to will be discussed at 
the end of this chapter, because it has non-trivial consequences. The real subtracted cross 
section is finite in d = 4 dimensions, as ti is allowed to approach infinity. 

Integrating the virtual operator V(r) without an infrared cutoff is more complex, con¬ 
sidering that there is an explicit integration over the splitting variables. Hence, we have to 
integrate this part in d = 4 — 2e dimensions analytically to extract the 1/e 2 and 1/e poles. 
The integrated virtual operator takes the form 


POO P 

/ dr V(r) = Y, dTl ~ to) = !(to) + K(f 0 ) , 

J to i J 


(3.8) 


where dT[ is the phase space integration of the additional parton. The decomposition of the 
integrated V(r) into two operators I (to) and K(to) is arbitrary. However, we choose I (to) 
to match the divergencies of the virtual amplitude, as it is customary. We emphasize this 
structure to illustrate that the parton shower naturally incorporates a subtraction scheme 
similar to the Catani-Seymour framework [14, 15]. In the case of initial state partons it 
is necessary to include additional collinear counterterms, denoted by P, needed for the 
renormalization of the parton distribution functions. Thus, the matched cross section reads 


- [F] = r [($,j p 0) ) + ( < i, m | p a)) + ($ m |[i( to ) + K(to) + p]|p^ ) ) 


ml 


+ /|S^ (jP|c/( ' F ’ io)l<I>m+l) 


(<wi/<>£U) - E(*-+iisrfWi - to) 


(3.9) 


For future reference, let us define the shorthands, as in Ref. [26]: 

($m|S) = (^ml^) + ($m\pW) + (<& m | [I(to) + K(t 0 ) + P] \p^ ) , (3.10) 

(<WI H) = (Zm+llp^+l) - Xl^m+llSzIp^)©^ - to) • (3.11) 

l 

The total cross section is then given by 

d NL °[l} = ^ J [d4> m ](l|4> m )(ch m |5) + J ld$ m+1 }(M$m + l)($m + l\H) , (3.12) 
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whereas including parton shower evolution amounts to the integrals 


ct[F] PS = — y f [d<S> m ](F\U(t F ,t 0 )\<S> m )($ m \S) 

m ■ J (3.13) 

+ , / [d$ m+ i)(F\U(t F ,t 0 )\$ m+ 1 )(<!> m+ 1 \H) . 

(m + 1)! J 

Matching is a two step procedure, which consists of first generating the samples according 
to Eq. (3.10) and (3.11), followed by the application of U(tp,to). 


3.3 Matching in the presence of singularities in the born approximation 

For processes with massless partons at leading order, the matching prescription as described 
in Eq. (3.9) must be slightly modified by the inclusion of generation cuts, as discussed for 
example in Refs. [58, 59]. A naive modification to the matching prescription would be to 
make the following replacements 

($m| S) ($ ro |S)lM{p, f}m) , (3.14) 

($m+i| H) ($ m+ 1 \H)F I ({p,f} m+1 ) , (3.15) 


where Fj(\p, f}\) is a jet function applied during the generation of events, on the momenta 
and flavours of 4 *a- Applying the parton shower to these ensembles shows that double 
counting is not removed. Indeed, substituting Eq. (3.14) and Eq. (3.15) in Eq. (3.13) we 
obtain 

p[F] PS = — [ [d«h m ](F|[/(t F ,t 0 )|4> m )(4> m |5)F / ({p,/} m ) 

m ' J /o -1 a\ 

i r (3.16) 

+ , -TTT / [d^m+l](F\U(t F ,to)\^m+l)(^m+l\H)F I ({p, f} m +l) • 

(m+ 1)! J 

Expanding the evolution operator, as given by Eq. (3.4), yields 


° [F]PS * / ld$mm$ m )($m\ [|^) + \p%) +P\p<®)] FdipJU) 


+ 

+ 


y j [d$ m+1 ](E|$ m+1 )($ m+1 |pS ) +1 )i 7 ’/({p,/Wi) 

dr (F|$ m+1 )($ m+1 |ft 7 (r)|$ m ) 


(m + 1) 

f ld$m} [d^m+l] f tF 

J ml (m + 1)! J to 

x ($m\p& ) )\F I ({pJ} m )-F I ({p,f} m+1 )\ +0(a 2 s ) , 


(3.17) 


where the I(to) + K(to) contribution of ( < 3? m |S') has been cancelled by the linear expan¬ 
sion of the Sudakov form factor. The offending term is present in the 3rd and 4th lines 
of Eq. (3.17). It does not vanish because Fj({p. f} rn ) is not equal to Fj({p, f} m +i) for 
non-singular configurations, despite the fact that their momenta and flavours are related 
according to {p, f} m + 1 = Ri({p, f}m , T/, xi)- The mismatch can be cured by enforcing the 
subtraction terms to fulfill Fj({p, f} m ), i.e. by modifying the real subtracted cross section 
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according to 


($m+l| H) 


($m+l| H) = - ^(Qm+llSllpWmti ~ i 0 )F / (Q I ({P, /}m+l)) , (3-18) 

l 

where we make use of the inverse momentum mapping Qi, as dehned in Eq. (2.30). Thus, 


Fi(Qi({p, f}m+ 1)) = Fjdp, f} m ) . (3.19) 


This modification allows us to compute cross sections with massless partons, at the same 
time introducing restrictions on the functional form of Fj. Using Eq. (3.18) in addition 
to Eq. (3.15) in the definition of the cross section in Eq. (3.13) and expanding the shower 
evolution yields 


a[F\ 


PS 


m\ 


+ 

+ 


[d<S> m ](F |<f> m )(<f> m | [| pW) + |p«) + P|p®)] FjdpJU) 
j [d$ m+1 ](F|$ m+1 )(# m+1 | /0 ^l)i ? /(fe/}m+l) 

dr (F|$ m+1 )($ m+1 |ft / (r)|$ m ) 


(m + 1)! 

f [d®m] [d®rn+ 1 ] f tF 

J ml (m + 1)! J to 

X (QmlpW) [l - F/({p, f}m+ i)l FjHp, f} m ) +0(a 


(3.20) 


The double counting is removed from Eq. (3.20), if the following condition is satisfied 


1 - F T ({p, f} m + 1) F({p, f} m+ i) = 0 , 


(3.21) 


where we have used the fact that (_F|3> m -|_i) ~ F({p, f} m+ 1 ). This is achieved if 


Fi({p, f} m + 1 ) = 1 , for F({p, f} m + 1 ) + 0 , (3.22) 

and can be understood as the generation cuts Fj(\p. /} m +i) being more inclusive than cuts 
that are applied on the final observable F({p, /} m +i). 

Our previous discussion of double counting relied on the fact that the splitting func¬ 
tions were the same in the NLO subtraction scheme and in the parton shower. However, in 
the case of initial state splittings, this is not the case due to mass effects and the presence 
of different PDFs. The parton shower PDFs are evolved as explained in Section 2.8 using 
modified leading order splitting kernels. On the other hand, the fixed order calculation 
employs NLO PDFs. We point out that the implied mismatch is of order O(a^) as long as 
both PDF sets are equal at some scale. This means that we use an NLO PDF parameteriza¬ 
tion as an input for the parton shower PDF evolution. A second source of mismatch is the 
treatment of mass effects in the initial state. In the case of the parton shower, the splitting 
kernels are modified as explained in Section 2.8, whereas the specifics of the treatment of 
mass thresholds depends on the PDF collaboration [60]. 

Finally, let us note that a last source of mismatch is due to scales in the strong coupling 
constant in the subtraction terms and in the shower. However, just like in the case of PDF 
evolution, it is of higher order. 
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3.4 On-shell projection 

As mentioned before, the Nagy-Soper parton shower allows for massive initial state partons 
and treats the charm and bottom quarks as massive. In an NLO calculation the charm 
quark is usually assumed massless, whereas the bottom quark is only present in the initial 
state if it is massless as well, which corresponds to the 5-flavour scheme. The practical 
application of the matching procedure requires, therefore, an on-shell projection for charm 
and bottom quarks. We use different procedures in the initial and final state cases. 


Initial state quarks 

We proceed iteratively by making one quark massive at a time. For definiteness, we will 
modify p a which is assumed to have pl = 0, whereas the mass, mj, of p b is arbitrary. The 
total momentum is given by 

Q=Pa+Pb , (3.23) 


where 

2 

rat 

Pa = rjaPA , Pb = VbPB H- ~PA , 

PbS 

with r] a and p b the momentum fractions, Pa/b the hadron momenta and s 
hadronic center of mass energy squared. This can be written as 

f m? \ 22 

Q = Pa d- } PA + VbPB , Q = SVaVb + m b . 

V VbsJ 

We shall determine rj a , such that 


a / / if] 2 / 

Pa = VaPA + —PB , P a = m a ^ 0 . 
VaS 


(3.24) 
2 (pa-Pb) the 


(3.25) 


(3.26) 


The new total momentum Q = p a + Pb an d its invariant mass then read 


rnt 


VbS 


Q = \Va-\ -- ) PA + \Vb + IT 2 - ) PB 


m„ 


VaS 


2™2 


Ao A 2 2 

Q = sp b pa + __ . +mi + mi 


SVbVa 


The requirement, Q z = Q , yields 


Va = Va I 2 


ml \! ( m a - Q 2 ) 2 - 4 mlml 

7TW7 T 77W) 


C IQ‘ 


: 2Q 2 


(3.27) 


(3.28) 


To ensure total momentum conservation, the complete final state has to be boosted in the 
^-direction. The rapidity of the boost, w, is 


e 


CJ 


Q 2 + ml - m 2 a 
2 ml 



which for m b = 0 reduces to 


Q 2 

Q 2 -ml' 


(3.29) 


(3.30) 
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Final state quarks 


In the case of a final state massless quark with momentum pi, we exploit the momentum 
mapping procedure of the Nagy-Soper parton shower to generate a non-zero mass mi. The 
new momentum is given by 

Pi = A Pi + 1 0 A + V Q , P 2 i = rn} , (3.31) 

zai 

where Q is the total final state momentum and 

- Q 2 , _ 

0,1 2pr Q ' 1 2pr Q ’ 

The parameter y is determined by requiring that the invariant mass of the total recoiling 
momentum, K = Q — pi, be preserved, i.e. K 2 = K 2 , with K = Q — pi. The result is 



V = k . (3.33) 

As a consequence, K and K are connected by a boost [16] given by 

(.K + ky(K + ky k»K v 

B IJ,/ = g> lu - 2---— + 2—— . 3.34 

(K + K) 2 K 2 ^ ’ 

The parameter A is determined by the on-shell condition p 2 = rn 2 . which yields 


A = y/(l + y) 2 - 4 a t bi 


(3.35) 


which is well-defined if 

Pi-Q > m t (^/Q 2 - . (3.36) 

In the center-of-mass frame of Q, this becomes 

E ‘ >ro, ( 1 “^*)- (337) 


Clearly, arbitrarily soft quarks cannot be projected on the mass-shell. 

In case of 2 —>• 2 processes, eg = 1 and the boost in Eq. (3.34) is singular. In that case 
one can use an alternative boost [16] 


with 


B ,,u = g^ + 


( K-n 
\K ■ n 


1 ^ n^fi v + ~ 1 ^ nV 


(3.38) 


n = 



n = 



(' Q ~ Pi) ■ 


(3.39) 
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3.5 Initial conditions for the parton shower 

As described in Section 2.7, the shower evolution emissions are strongly ordered in the 
parameter A 2 , which is related to the shower time t as follows, 


A? 


P -t - tlL A 2 - 

Ql' 


±Pm+ 1) 2 ~ mf 

2pi ■ Qo 


Qo > 


(3.40) 


where pi is the emitter momentum and Qo is the total final state momentum before emission. 
The emitter momentum after emission is pi, while p m +1 is the momentum of the emitted 
particle. The parton shower evolution starts at an initial time to, which has to be determined 
by the user. Since the choice of to has a noticable influence on the differential distributions, 
it is necessary to provide a sensible prescription, which is the purpose of this Section. 

Consider first a process at leading order, e.g. pp —>• tt. The parton shower will generate 
additional radiation. However, the resulting hard jets are described very poorly, because 
the parton shower description is only valid close to infrared limits. As a consequence, the px 
distribution of the top quark pair is not reliably reproduced for large values of the transverse 
momentum. Since the initial time to restricts the available phase space for parton emissions, 
it will have a strong influence on the distribution, as reported, for example, in Ref. [61]. 
Notice that the situation may be improved with merging of different multiplicity samples, 
see Ref. [62]. 

Matching at NLO has the advantage of providing one hard emission exactly through the 
matrix element. According to our modified quantum density matrix, we removed all parton 
shower contributions to the first emission. Therefore, corrections to the px distribution are 
formally of higher order, even though they can be a priori large [57, 63]. The ultimate goal 
is to choose to such that the description of the high px regime is close to that of the NLO 
fixed order approximation. This can be achieved by restricting the parton shower phase 
space to allow only splittings with virtualities lower than the splittings already present in 
the hard matrix element. 

We choose to to approximate the virtuality of the hard process on an event-by-event 
basis. Inspired by Refs. [27, 64], we define 


e to = min 

*7O' 


2Pi • Pj } 

Qo J 


(3.41) 


where pt and pj are pairwise different external momenta of the process considered. 

For the real radiation contribution, we use the subtraction terms Eq. (3.18) to determine 
to as follows 


1. For each subtraction term 


(4? m+1 |S J |pW)0(tj-t o )F / (Q i ({P,/}m + i)) , (3.42) 

determine its individual to according to Eq. (3.41) while ti according to Eq. (3.40). 
Retain ti if ti > t 0 . 
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2. If there is at least one t[, then evaluate to according to 

e~ to = min (e _ti ) , (3.43) 

where the minimum runs over all ti found in the previous step. 

3. Otherwise, apply Eq. (3.41) to the [m + l)-particle kinematics. 

Let us note that the separate determination of to hi step 1 for each subtraction term, is 
necessary to avoid double counting. 

3.6 A summary of ambiguities 

In this Section, we recapitulate by listing the intrinsic uncertainties introduced by the 
Mc@NLO matching formalism in combination with the Nagy-Soper parton shower. 

Parton distribution functions: Parton distribution functions are evolved differently 
in the NLO calculation and in the parton shower. Nevertheless, we pointed out in 
Section 3.3, that the evolution itself is of higher order. Thus, the NLO accuracy is 
obtained as long as the evolutions share a common point, e.g. at the low scale. 

Parton masses: A generic feature of the Nagy-Soper parton shower is that it requires 
massive bottom and charm quarks in the initial state. However, the fixed order 
calculation involves massless charm quarks at least. Although the implied mismatch 
is power suppressed, it is necessary to introduce masses for the relevant quarks in the 
fixed order sample used as input for the shower. In Section 3.4, we have provided a 
possible algorithm. Nonetheless, other choices might influence the final result. 

Initial state shower time: The choice of to in the parton evolution is, to a large extent, 
arbitrary. The only requirement is that the NLO prediction be recovered for hard 
emissions. This can be achieved by several different choices of to, with either fixed 
or configuration dependent values. In Section 3.5, we have presented our choice, but 
others are possible. Nevertheless, it is recommended to vary to by rescaling it by a 
common value within some reasonable range, see Section 5.1. 

4 Implementation 

In this Section, we present a realization of the matching scheme within the HELAC-NLO 
framework [65]. Due to the restricted functionality of the Nagy-Soper parton shower present 
in DEDUCTOR version 1.0.0 [24], which we used as basis for practical studies, our implemen¬ 
tation is simplified as far as spin and colour treatments are concerned. Since DEDUCTOR 
uses spin averaged splitting functions, we only provide unpolarized event samples for show¬ 
ering. On the other hand, we supply only leading colour events even though DEDUCTOR 
works with the LC+ approximation and may shower non-diagonal colour configurations. 
While proper colour matching does not present conceptual challenges, it requires some ad¬ 
ditional programming effort. The latter will be necessary once DEDUCTOR has full colour 
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functionality. We believe that, for a first study, our simplification is justified. More details 
on our approach can be found in Appendix A. Finally, we note that, at present, Deductor 
lacks unstable particle decays and a hadronization model. 

4.1 Modifications in Helac-Dipoles 

The implementation of the matching scheme in the Helac-Nlo multi-purpose event gen¬ 
erator only concerns HELAC-DlPOLES [33, 35]. We base our work on the previous imple¬ 
mentation [33], where the majority of the parton shower operators Tij and V have already 
been included. 


Momentum mapping for initial state splitting 

The momentum mapping implemented in HELAC-DlPOLES is based on Ref. [16]. For rea¬ 
sons, which we shortly discussed in Section 2.4, the mapping implemented in DEDUCTOR 
is slightly different [19]. Here, we reproduce the necessary formulae used in the new version 
of Helac-Dipoles. 

We start from a set of momenta {p} m and consider the splitting p a —>• p a + p m+ The 
incoming momenta are given by 


Pa = fjaPA , Pb = VbPB , 


(4.1) 


where Pa and Pb are the hadron momenta. After the splitting, there is 

1 

Pa = ~Pa , Pb=Pb , 

z 


and the momentum p m + 1 °f the emitted particle is parametrized by 


(4.2) 


Pm +1 = XaPa + %bPb + , (4.3) 

where 

|/c_l| 2 = 2 x a x b (p a -p b ) , x a = - - 1 - y , x b = zy . (4.4) 

The phase space variables y, z and (f> are defined as 

r _ rja _ sfjgfjb _ Q 2 

Pa SpaPb (Pa + PbY ’ 

_ (p a -pm+l) 2 _ 2 p a • p m+l (4.5) 

V 2(Pa ■ Pb) Q 2 

<t> e [0, 2 t r] , 


where Q = p a + Pb ~ Pm+l = Pa + Pb- T° ensure momentum conservation the remaining 
final state particles have to be boosted, pi = Bpi, with the Lorentz transformation [22] 


B^(u,v ± )=g T + ^K+Cm 


Pa'Pb 


+ 


Pa 'Pb 


~Pb v ±\ ~ 


P~ V 


2 PbPb 


‘ Pa-Pb 


(4.6) 
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where 


uu U u PtPb+PbPa 
y_L — y 


Pa -Pb 


1 


e =- x a = l + y , 

z 


V± = - 


V 2 Pa ■ Pb 


kj_ . 


(4.7) 


In the case, where pb is the emitter one simply exchanges p a and pb in the equations above. 
This momentum mapping is implemented in the I and KP operators of the HELAC-DIPOLES 
package, while the modification only concerns initial-final state interference contributions. 
For the real subtracted cross section, the inverse transformation has to be applied 


Pa = zp a , Pb=Pb , 

Pi = B~ 1 p i (i = 1 , ...m) , 


(4.8) 


with 


(B-T = ar + + ^ m 

Pa'Pb 


+ 




(4.9) 


Pa'Pb Pa 'Pb 


Dynamical phase space restriction 

A second modification to the existing subtraction scheme is the implementation of a dynam¬ 
ical cutoff @(t; — to) on the available dipole phase space, which represents the ordering of the 
emissions in the parton shower evolution. The new subtraction terms read (see Section 3.2) 


POO _ 

/ dr Ui(t) = SjQ(ti - t 0 ) , 
Jto 1 


(4.10) 


where the splitting operator S i is already available in HELAC-DlPOLES. The shower time 
ti is given by (see Section 3.5) 

_ t \(pi ± p m+ i) 2 - mf 

6 1 — - 

2 Pi ■ Q o 

and can be reconstructed from the {p, /} m +i kinematics. The numerator is given by real 
radiation momenta, whereas the denominator 2 pi ■ Q o is given by the momenta before the 
splitting. However, the denominator can also be reconstructed by making explicit use of 
the momentum mappings. 

This phase space restriction is also implemented in the integrated dipoles. In order to 
use as much existing code as possible, we implement it in the form of 


(4.11) 










Therefore, we can calculate the finite remainder in d = 4 dimensions and subtract it from 
the complete expression. We follow closely the outlined semi-numerical strategy presented 
in [33], 

We note that, even though the restriction is necessary for matching, it can be also used 
in fixed order NLO calculation as an effective cutoff on the subtraction phase space similar 
to ct max in the Catani-Seymour subtraction scheme [66-68]. 


4.2 Monte Carlo techniques 


With the modifications described above, it is possible to use Helac-Nlo to generate an 
event sample ready for showering with DEDUCTOR. We produce events subprocess by 
subprocess. First, we use HELAC-1LOOP [69] to obtain a set of unweighted leading order 
events with the virtual contributions 


^iii.Pi /}m) 


1 + 


j{P,f}m\Pm) 
({ P,f}m\Pm ) 


where only the finite part of Ip'm) is included. At this point 


(4.13) 


cr[LO + V] 


cr[ LO] 
N 


N 

J2 w i ■ 

i= 1 


(4.14) 


This set of weighted events is subsequently reweighted using HELAC-DlPOLES in order to 
include the parton shower virtual operator. This corresponds to taking into account the 
integrated subtraction terms. The weights become 


^iiiPi /}m) 


1 + {{P,f}m\Pm) + ({p, + K(£q) + P|pl?) 

i{Pif}m\Pm) ({ P,f}m\Pm ) 


(4.15) 


yielding 

a[LO + V + I + KP] = . (4.16) 

i— 1 

Notice that the integrated subtraction terms are Monte Carlo integrals on the phase space 
dTi of the additional unresolved parton. In order to obtain a good approximate of this 
integral, we sample several points for a fixed born phase space point {p, f} m ■ This proved 
to be advantageous in case of large cancellations in V(r). 

The real radiation events are generated separately with HELAC-DlPOLES, which was 
extended to provide unweighted events with positive and negative weights, ±1 [70]. Notice 
that unweighting is possible, because both the real radiation and the respective subtraction 
weights correspond to the same phase space point as described in Section 3.2. For each 
accepted event, we choose the most probable diagonal colour flow with the colour weight 


CwiAPi /) c }m+l) 


{{pj, cU+iIpS^) 

£{c'} m+1 ({p> /> c 'Wl |Pm+l) 


(4.17) 


The generated events are stored in a Les Houches hie [71], with the initial parton shower 
time to assigned to the variable SCALEUP. 
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4.3 Interface to Deductor 


Events generated by Helac-Nlo are transferred to DEDUCTOR. This requires an on- 
shell projection for charm and bottom quarks (see Section 3.4) and the determination of a 
starting colour configuration for each event. 

On-shell projection 

The on-shell projection of Section 3.4 is applied iteratively to each charm and bottom 
quark. In the case of an initial state transformation for two massless quarks bound to 
become massive, the order of longitudinal boost is relevant. Indeed, the first boost has a 
lower rapidity than the second as can be proven by inspection of Eqs. (3.29) and (3.30). We 
choose the order at random in order to reduce the systematics. In case the projection fails, 
which is only possible for soft bottom and charm quarks in the real radiation contribution, 
the event is rejected. Vetoing such emissions may only induce a negligible modification of 
the cross section. Otherwise, effects, which are not under control, would be substantial. 
This problem has to be studied case-by-case. 

Colour configurations 

A colour flow generated as described before in the Les Houches [71] format must be trans¬ 
lated to the internal representation of DEDUCTOR in terms of colour strings, see Section 2.1. 
Notice that only one colour flow is needed in the leading colour approximation per event. 
A colour flow is given by a list of two colour indices. Thus, the f-th particle carries the 
following pair of indices 


q(i) (colourlfi] , 0) 

q(i) —> (0 , colour2[i]) (4-18) 

g(i) —>• (colourl [i] , colour2[i]) . 

In DEDUCTOR, on the other hand, quarks, anti-quarks and gluons are represented by 

q(i) -»■ [Q , next [i]) 

q{i) (prev [i] , A] (4.19) 

g(i) -»■ (prev [i] , next[i]) , 

where Q denotes the beginning and A the end of an open colour string. prev[i] refers to 
the particle index in the event, which is to the left of particle with index i. In the same 
way next [i] refers to the right partner on the colour string. The algorithm to translate 
colour flows into colour string configurations is 

1. Swap the flavour and the colour indices for initial state partons. 

2. Enumerate all partons starting from —1, where —1 and 0 are reserved for the initial 
state partons. Colour-neutral particles are enumerated with numbers less than —1. 

3. Iterate over all particles with index i applying 
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(a) If (colourl [i] 7 ^ 0 and colour 2 [i] = 0) then prev[i] - Q. 

If (colourl [i]= 0 and colour2[i] 7 ^ 0) then next [i] A. 

(b) Iterate over all particles with index k 7 ^ i applying 

If (colourl [i] 7 ^ 0 and colourl [i] = colour 2 [k]) then next [i] —k; prev [k] - i. 


5 ttj production at the LHC with next-to-leading order matching 

In this Section, we present results for pp —>• ttj + X production at next-to-leading order, 
obtained with Helac-Nlo, matched with the Nagy-Soper parton shower as implemented 
in DEDUCTOR. The NLO QCD corrections to the considered process have been previously 
computed in Ref. [72-74]. Matching to a parton shower has been first considered in Ref. [75, 
76] using the POWHEG method [28, 30]. 


5.1 Setup 


The results for ttj production are presented for pp collisions at the LHC with a center- 
of-mass energy of 8 TeV. The top quark is assumed to be stable and its mass is set to 
nit = 173.5 GeV, while the charm and bottom quarks are considered to be massless at fixed 
order. Results are obtained using the Mstw2008nlo PDF set [77] with five active flavours 
and the corresponding two-loop running of the strong coupling. We set the renormalization 
and factorization scales to the top quark mass, pr = Pf = mt, and the starting shower 
time to 


e = min 

*7A? 


• Pj \ 

p 2 T Ql J 


(5.1) 


where pi and pj are external momenta, Q 0 is the total final state momentum and px = 1 for 
the central prediction as explained below. Partons with pseudorapidity \rj\ < 5 are clustered 
using the anti-fey jet algorithm [78], with the separation parameter R = 1. The resulting 
jets are sorted in decreasing order of px- We require the tagged jets to have transverse 
momentum of px > 50 GeV and rapidity in the range of |y| < 5. 

We restrict our analysis to the perturbative parton shower evolution. Decays of unstable 
particles, hadronization and multiple interactions are not taken into account. The parton 
shower treats the charm and bottom as massive particles with masses m c =1.4 GeV and 
mi = 4.75 GeV. We provide the Mstw2008nlo PDF set at pp = 1 GeV as the starting 
point for the evolution in DEDUCTOR. We also use the corresponding two-loop running of 
a s , and restrict the parton shower to the leading colour approximation. Therefore, results 
presented in this section are accurate up to 0(1/IV,?). 

In order to address the theoretical uncertainties, we investigate the scale dependence of 
cross sections and distributions on the unphysical scales pr, pf and the rescaling parameter 
px- Here, pr is varied simultaneously with pp between pr,f = Rv/2 and pr } f = 2m;. 
On the other hand, both the central value of px, Pto , and its variation range require a 
more thorough discussion, because this parameter is specific to the Nagy-Soper parton 
shower. Guidelines for a suitable choice of to have already been presented in Section 3.5. 
In practice, we consider exclusive distributions, which are especially sensitive to the parton 
shower effects, e.g. the transverse momentum, px, spectrum of the ttj\ system or of the 
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Figure 10. Differential cross section distributions as a function of the transverse momentum of 
the ttj i system (left panel) and of the second hardest jet (right panel) for pp —> ttj + X at the LHC 
with y/s = 8 TeV. Comparison between the NLO result obtained with Helac-Nlo and the results 
produced by matching NLO predictions to Deductor for three different choices of the starting time. 
The scale choice is pf = pr = mt ■ The lower panels display the relative deviation from the fixed 
order result. Notice, that to = 0 corresponds to unrestricted shower radiation without relation to 
the kinematics of the underlying event. This choice should not be used in practice and it is shown 
for illustration purposes only. 


second jet, j 2 , which are equivalent at NLO. For large pp values, the fixed order prediction 
is reliable and we would not like the shower to introduce substantial shape differences there. 
This can be achieved by a suitable choice of pp 0 as shown in Fig. 10, where the lower panel 
shows the relative deviation from the fixed order result defined as 

Aa _ (jNLO+ps _ ct nlo 

V = 

The pp spectra of ttj 1 and j '2 coincide at NLO because of momentum conservation. 
After including shower effects, they deviate since the ttji system recoils against several 
jets. In addition, the parton shower generates corrections to the structure of the second 
jet. We can see that for to = 0 and pp = 2 the parton shower overshoots the tail of 
the transverse momentum spectra. The large higher order corrections can be explained by 
the exponentiation of non-singular emissions, as already reported in Ref. [57] in the case 
of Higgs production. Decreasing the value of pp helps to recover the NLO predictions. 
We conclude that an appropriate choice of the central value of the starting time rescaling 
parameter is 

MTO = 1 . (5.3) 

We expect that the choice of pp will affect all exclusive distributions. We present two 
examples in Fig. 11: the rapidity of the second jet, y(j 2 ), and the number of hard jets. 
Finally, after fixing the central value of pp, we choose to vary it between pp = ppo/2 and 


(5.2) 
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Figure 11. Differential cross section distribution as a function of the rapidity of the second hardest 
jet (left panel) and inclusive jet cross sections (right panel) for pp —> ttj+X at the LHC with y/s = 8 
TeV. Comparison between the NLO result obtained with Helac-Nlo and the results produced by 
matching NLO predictions to Deductor for three different choices of the starting time. The scale 
choice is fiF = Hr = cnt- The lower panels display the relative deviation from the fixed order result. 


Ht = 2/iTO- This initial time variation is the dimensionless analogue of the resummation 
scale variation originally introduced in Refs. [79-83] and first used in the context of event 
generators in Ref. [84], 


5.2 Results 


The fixed order cross section at NLO obtained from the HELAC-NLO framework, including 
the theoretical error estimated via scale variation, reads 


NLO_ 

a pp^ttj+X 


86.04 


+ 5.10 ( + 6 %) 
- 11.41 (- 13 %) 


pb. 


(5.4) 


In order to obtain a reliable result including showering effects, we study the dependence of 
the cross section on the generation cut, denoted by p^. Table 1 contains values obtained 
at pp = hr = mt and pp = 1- We see that for all values of p^ but the last one, for which 
the generation cut is equal to the analysis cut, the total cross section is compatible with the 
fixed order prediction. For the following study, we choose a generation cut of pp(j l) > 30 
GeV. Varying the renormalization/factorization scale between mt /2 < pr^f < 2mt and 
the initial shower time rescaling parameter between 1/2 < pp < 2, the total cross section, 
together with its uncertainties, obtained after applying the parton shower is 


NLO+PS 

pp—yttj+X 


86.11 


+ 4.38 ( + 5 %) 
- 10.88 (- 13 %) 


[scales] 


+ 0.80 (+ 1 %) 
+ 2.17 (+ 3 %) 


[PS time] pb , 


(5.5) 


where the upper, (lower) values are given for /iR,F = mt/ 2 (2 mt) and pp = 1/2 (2). The 
scale dependence of the total cross section, taken very conservatively as a maximum of the 
upper and lower results, is 13% or 9% after symmetrization. The dependence on the initial 
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j/f jt [GeV] 

CS [P b ] 

e [%o] 

5 

86.51 ±0.21 

2.4 

10 

86.26 ±0.17 

2.0 

15 

86.22 ±0.14 

1.6 

30 

86.11 ±0.13 

1.5 

40 

86.01 ± 0.08 

0.9 

50 

84.58 ± 0.07 

0.8 


Table 1. Total cross section for pp —> ttj+X at the LHC with y/s = 8 TeV, together with statistical 
and relative errors, for different values of the generation cut. Results are produced by matching 
Helac-Nlo predictions to Deductor. The cross section is calculated for pp = pr = mt and 
Pt = 1 - 
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Figure 12. Differential cross section distributions as a function of the transverse momentum of the 
ttj i system (left panel) and of the hardest jet (right panel) for pp —> ttj + X at the LHC with y/s = 8 
TeV. Results are produced by matching Helac-Nlo predictions to Deductor. The uncertainty 
bands depict scale and initial shower time variation. The lower panels display the corresponding 
relative deviation from the central value, separately for pr,f a,nd pp. 


shower time is, by comparison, negligible. The situation is quite different when it comes 
to differential distributions as can be observed in Figs. 12 and 13, where the transverse 
momentum of the ttj\ system, the first and second hard jets together with inclusive jet 
cross sections are given. The variation bands for Pr,f an d pp in Figs. 12 and 13 have been 
obtained using the following sets of three parameter values: Pr,f = {r7Tt/2, mt, 2mt} and 
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Figure 13. Differential cross section distribution as a function of the transverse momentum of 
the second jet (left panel) and inclusive jet cross sections (right panel) for pp —> ttj + X at the 
LHC with y/s = 8 TeV. Results are produced by matching Helac-Nlo predictions to Deductor. 
The uncertainty bands depict scale and initial shower time variation. The lower panels display the 
corresponding relative deviation from the central value, separately for pr^f and pp. 

pp = {1/2, 1, 2}, respectively. The lower panels of Figs. 12 and 13 display corresponding 
relative deviations from the central value, separately for p^F and pp. 

We start the discussion with the transverse momentum of the pp(ttji) system, which is 
presented in Fig. 12 (left panel). At leading order, this observable is zero due to momentum 
conservation. When real emission contributions at the NLO level are included, this observ¬ 
able diverges as the transverse momentum of the entire system goes to zero. Therefore, it 
can only be reliably described by the fixed order calculation in the high pp region. However, 
including the parton shower, the low pp behavior is altered strongly by the Sudakov form 
factor as can be seen in Fig. 12. Indeed, for low values of the transverse momentum, the 
distribution is generated mostly by the parton shower. The reason is that the real radia¬ 
tion contribution, which is responsible for the divergent behaviour at fixed order receives 
subtractions, which match the singular behaviour for pp —>• 0. These subtractions belong 
to the same bin as the real radiation events themselves, contrary to the fixed order calcu¬ 
lation, where they belong to the zero bin. Thus, the subtracted real radiation sample has 
low weight contributions for low pp. In consequence, we observe a moderate dependence 
on pp in this region, which reaches factors of 1.5 at the lower end of the spectrum. This 
dependence decreases down to just a few percent around 30 GeV, whereas for moderate 
values of pp(ttji) it is at the level of 20% — 45%. The presence of a minimum dependence 
on pp is due to a crossing of the distribution for pp = 1/2, which dominates at low pp 
with the distribution for pp = 2, which dominates at high pp. The reason for a larger cross 
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section for low values of px and px is that the parton shower generates low px radiation 
only barely shifting events with zero px to non-zero values. On the other hand, high fix 
results in radiation of high px partons, which shift events to high px values. The situation 
is reversed when it comes to the renormalization and factorization scale dependence. Here, 
visible deviations from the central value occur once the matrix element is present. They 
grow substantially up to almost 80% at the end of the spectrum. This can be explained 
by the fact that the variation of Pr,f is only implemented in the matrix element, while the 
shower does not depend directly on those scales but rather on to- 

The px(ji) distribution of the hardest jet, which is given in Fig. 12 (right panel) shows 
a rather constant and small dependence on both parameters, Pr/f an d pxi as it has NLO 
accuracy. The Px{j 2 ) distribution of the second jet, which is presented in Fig. 13 (left 
panel), also shows a rather constant scale dependence with a somewhat larger variation 
range, as it is only LO accurate. 

Finally, inclusive jet cross sections are shown in Fig. 13 (right panel). As expected, the 
NLO cross section with exactly one jet, which is given in the first bin, is rather insensitive 
to px- Its theoretical error is at the 12% level. The px dependence is slightly larger in the 
second bin, where the two jet cross section, correct only at the LO level, is stored. Also 
here, the theoretical error increases up to 16%. Starting from the third bin, cross sections 
are described via the shower evolution alone, therefore, fairly large variations can be noticed 
for both parameters, px and pR t F- For example, the scale dependence for the cross section 
with five jets is found to be around 35%. In addition, cross sections for higher multiplicities 
are severely affected by the cutoff on the parton shower emission phase space. This strong 
dependence can be partially cured when matrix elements for higher multiplicities are merged 
together and matched to the shower evolution as described e.g. in Refs. [62, 85] for the LO 
and in Refs. [64, 86] for the NLO case. 

5.3 Comparison with other Monte Carlo event generators 

In this section we compare our implementation of the NLO matching to the Nagy-Soper 
parton shower from DEDUCTOR, as implemented in HELAC-DlPOLES, with other matching 
procedures and shower programs. To be more specific, we use: 

1. the POWHEG-BOX [30] implementation of pp —>• ttj + X [76] in conjunction with the 
PYTHIA 8.1 ]45] Monte Carlo program version 8.183, with the transverse-momentum 
ordered shower (dubbed POWHEG+PYTHIA8), 

2. the automatic solution of aMc@NLO [59] including Pythia 6.4 [44] (version 6.428) 
with the virtuality-ordered or mass-ordered shower (dubbed aMc@NLO+PYTHIA6Q), 

3. the automatic solution of aMc@NLO together with PYTHIA 8.1, once more with the 
transverse-momentum ordered shower (dubbed aMc@NLO+PYTHIA8). 

Once again, we do not include top quark decays, hadronization and multiple interactions. 
Thus, we are only comparing the perturbative evolution of different parton showers. In 
each case the default setup of the programs is used. 
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The total cross sections together with their theoretical errors are 

a pp°fff+x (aMc® Nlo-| Pythia 6 Q) = 84.85[scales] pb , (5.6) 

C-?i P +x( aMc@NLO + PYTHIA8 ) = 89 - 55 -i 5 .’ 4 i (— 17 %) [scales] pb , (5.7) 

(Powheg+Pythia8) = 89.12+2' 22 ( +f 0 g [scales] pb . (5.8) 

We observe that all three calculations, although based on different shower ordering variables, 
give compatible results and agree within 4% with the HELAC-NLO+DEDUCTOR result from 
Eq. (5.5). We note that, in case of PYTHIA8, the central prediction for our process is larger. 
Moreover, the scale dependence after symmetrization is below 13% for all cases but the 
POWHEG+PYTHIA8, where it is slightly larger i.e. of the order of 20%. 

In the next step, we extend our comparison to differential distributions and start with 
observables that are rather insensitive to parton shower effects. In Fig. 14, the transverse 
momentum and rapidity distributions of the top quark and the first jet are presented. 
For each observable the NLO result obtained with HELAC-NLO is plotted together with 
results produced by matching various NLO predictions to different parton showers. The 
lower panels display the relative deviation from the next-to-leading order result. All parton 
showers reproduce the corresponding NLO result correctly, and slight deviations from the 
fixed order calculation can only be seen in the tails of the distributions due to smaller 
statistics. 

The same conclusions can be drawn from Fig. 15, where we show the invariant mass 
of the tt pair and the angular separation in the rapidity-azimuthal angle plane between the 
top and the anti-top quark 

ARij = sj(yi - yj ) 2 + (fa - (pj) 2 . (5.9) 

Overall, we observe that results obtained within the HELAC-NLO+DEDUCTOR framework 
are consistent with those of the other Monte-Carlo event generators used in this study. 
Let us emphasize, that not only various showers but also different matching schemes with 
distinct systematic uncertainties are examined here. 

To assess potential differences among showers and matching procedures used in our 
analysis, we turn to observables, which are sensitive to the initial conditions of the parton 
evolution. As an example, we present in Fig. 16 the differential cross section distribution 
as a function of the transverse momentum of the ttji system. For clarity we give both a 
linear-scale (left panel) and a log-scale (right panel) version of the plot. We observe that 
the matching with PYTHIA8 yields similar results for both Mc@Nlo and POWHEG, i.e. 
the hardest spectra and a lack of agreement with the fixed order prediction in the region 
where hard, well separated partons are produced. We have found discrepances even up 
to 120%. This suggests, that for the process at hand, initial condition of the transverse- 
momentum ordered shower from PYTHIA8 should be adjusted to restrict further the shower 
phase space where emissions become hard. On the other hand, aMc@NLO+PYTHIA6Q 
nicely reproduces the NLO differential cross section in the high pr region. We also note 
that it produces a softer spectrum. As for the HELAC-NLO+DEDUCTOR case, as already 
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Figure 14. Differential cross section distributions as a function of the transverse momentum of 
the top quark (left upper panel) and of the first jet (left lower panel) as well as the rapidity of 
the top quark (right upper panel) and of the first jet (right lower panel) for pp —> ttj + X at the 
LHC with ^/~s = 8 TeV. Comparison between the NLO result obtained with Helac-Nlo and results 
produced by matching various NLO predictions to different parton showers. The scale choice is 
fip = pp. = ffh- The lower panels display the relative deviation from the fixed order result. 


explained, the parameter to was chosen to preserve the NLO shape in the tail of the pp 
spectrum of the ttji system. 

In Fig. 17, the kinematics of the second jet is analysed. To be more precise, distri¬ 
butions in the transverse momentum and in rapidity are given. There is a good agree¬ 
ment between the spectra of aMc@NLO+PYTHIA6Q and Helac-Nlo+Deductor, with 
a somewhat narrower rapidity spectrum in the case of HELAC-NLO+DEDUCTOR. For 
aMc@NLO+PYTHIA8 and POWHEG-I-PYTHIA8, a difference in normalization can be ob¬ 
served, but the differential K -factor remains flat. This can be explained by larger emission 
rates that increase the amount of radiation. 
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Figure 15. Differential cross section distributions as a function of the invariant mass of the tt 
pair (left panel) and of the NR t i (right panel) for pp —> ttj + X at the LHC with i /s = 8 TeV. 
Comparison between the NLO result obtained with Helac-Nlo and results produced by matching 
various NLO predictions to different parton showers. The scale choice is pf = Pr = m t ■ The lower 
panels display the relative deviation from the fixed order result. 




P T (ttii) [GeV] log 10 p T (ttj-|) 


Figure 16. Differential cross section distributions as a function of the transverse momentum of 
the ttji system for pp —> ttj + X at the LHC with y/s = 8 TeV. Comparison between the NLO result 
obtained with Helac-Nlo and results produced by matching various NLO predictions to different 
parton showers. The scale choice is pf = /i# = m*. The lower panels display the relative deviation 
from the fixed order result. 

In Fig. 18, the differential cross section as a function of the transverse momentum of 
the tt pair is presented. We observe a good agreement among all programs for pp > 50 
GeV i.e. above the analysis cut. Below this value the prediction is only leading order 
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Figure 17. Differential cross section distributions as a function of the the transverse momentum 
(left panel) and rapidity of the second jet (left panel) forpp —> ttj + X at the LHC with y/s = 8 TeV. 
Comparison between the NLO result obtained with Helac-Nlo and results produced by matching 
various NLO predictions to different parton showers. The scale choice is hf = hr = W(. The lower 
panels display the relative deviation from the fixed order result. 
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Figure 18. Differential cross section distribution as a function of the transverse momentum of the 
tt pair (left panel) and inclusive jet cross sections (right panel) for pp —> ttj + X at the LHC with 
y/s = 8 TeV. Comparison between the NLO result obtained with Helac-Nlo and results produced by 
matching various NLO predictions to different parton showers. The scale choice is /if = h-R = m*. 
The lower panels display the relative deviation from the fixed order result. 

accurate and thus strongly depends on the initial shower conditions. Also shown in Fig. 18 
are inclusive jet cross sections. In the first (second) bin the exclusive cross section is 
accurate at next-to-leading order (leading order). Starting from the third bin on, cross 
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Framework 

jet] 

a [>2 jets] 

,j[>3 jets] 

a [> 4 jets] 

ct [>5 jets] 

Helac-Nlo+Deductor 

86110 

17204 

2780 

362 

38 

aMC@NLO+PYTHIA6Q 

84850 

15030 

1230 

106 

7 

aMc@NLO+PYTHlA8 

89556 

21872 

3377 

372 

32 

Powheg+Pythia8 

89121 

23744 

5458 

992 

129 


Table 2. Cross sections (in fb) for the inclusive jet rates at the LHC with y/s = 8 TeV, according 
to the default settings of the various codes. 


Framework 

C| 

IV 

IV 

O'!- 3 ] /crl— 2 1 


crl— 5 ! /crl— 4 ] 

Helac-Nlo+Deductor 

0.20 

0.16 

0.13 

0.10 

aMC@NLO+PYTHIA6Q 

0.18 

0.08 

0.09 

0.07 

aMc@NLO+PYTHlA8 

0.24 

0.15 

0.11 

0.09 

Powheg+Pythia8 

0.27 

0.23 

0.18 

0.13 


Table 3. Cross section ratios for (n + 1 )/n inclusive jet rates at the LHC with y/s = 8 TeV, 
according to the default settings of the various codes. 


sections are only described with leading-logarithmic accuracy via the parton shower alone. 
Results for exclusive cross sections with JVj e t s = 1 and A(j ets = 2 are in agreement with our 
observations for other observables. The picture is vastly different for all showers, however, 
when iVjets > 3. In this case, the Powheg+Pythia 8 framework (aMc@NLO+PYTHlA6Q) 
produces the highest (lowest) number of hard jets. 

Finally, cross sections (in fb) for inclusive a(tt + n jets) rates are presented in Table 2. 
Table 3 contains cross section ratios, i.e. a(tt + n + 1 jets)/<r(ff + n jets). 

6 Conclusions 

In this publication, we have presented a next-to-leading order matching scheme for the 
Nagy-Soper parton shower. We based our construction on the original Mc@Nlo approach. 
Besides the general formulation, we have performed real simulations for top-quark pair pro¬ 
duction in association with a jet at the LHC, using an implementation within the framework 
of the public codes Helac-Nlo and DEDUCTOR. 

Our general conclusion is that the combination HELAC-NLO+DEDUCTOR is able to 
provide results for non-trivial processes, which remain in reasonable agreement with other 
Monte Carlo systems. Indeed, for observables, which are rather insensitive to showering 
effects, the differential cross sections are in very good agreement between different programs 
and the fixed order NLO calculation. Inevitable differences for infrared sensitive observables, 
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on the other hand, seem to be justifiable in size in the sense that all predictions have 
overlapping uncertainty bands. Of course, further studies are needed here. 

We note that our simulations are at the same level of logarithmic precision as those of 
others. For now, we are only correct as far as the leading behaviour is concerned. Future 
developments in DEDUCTOR will allow us to include soft-gluon intereference effects yielding 
next-to-leading logarithmic accuracy. This step still requires some improvements of our 
implementation in HELAC-DlPOLES. In particular, it will be necessary to transfer colour- 
configuration information exactly, as opposed to the current leading-colour approximation. 
We leave this to future work. 

Finally, we expect that the Nagy-Soper parton shower will open new opportunities for 
understanding parton shower systematics for processes with non-trivial colour exchange. 
Here, the road is still long, as we must remember that, ultimately, hadronization models 
must be included. The latter, however, require tuning to the shower. There are also 
interesting problems in merging different multiplicity samples generated by our software. 
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A Simplified matching scheme 

In this Appendix, we discuss the simplified matching prescription that we have implemented 
starting from spin and colour averaged amplitudes. We also prove that it yields NLO 
accuracy at leading colour, as long as the observable is not sensitive to spin correlations. 
The simplifications match the current DEDUCTOR functionality and relate to a spin and 
colour averaged quantum density 

\p) = X! “[ f f}m]\{p, f}m)({p, f}m\p ) , (A.l) 

m ' J 

where 


{{P,f}m\p) = y. e ({s, c} m |{s', c'} m ) p({p, /, s', c', s, c} m ) . (A.2) 

R}m,RTm lUm rn 

Averaging over spin and colour reduces tremendously the complexity of the calculation. 
However, the parton shower requires a colour configuration for each phase space point 
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fa, f}m in order to perform the evolution. Starting from colour averaged weights (fa, f} m \p) 
one can recover leading colour correlations, (fa,/, c} m |p), by including a colour weight 
C w ({p, defined by 


fau (fa) f 1 c}m) 


(fa) /) c} m \Pm) 
E{c} m (fa /> C}m| Pm) 


The simplified matching prescription accounting for generation cuts is then 


(A-3) 


a[F \ = I \d{p, f, c} m \ /? c } m )Cfafa, /, cfa)(fa, /} m |fapW)ifafa, f}m) 

J ml 

+ [ [d{P ' f, " } Z +l] ( F \U(t F , to)\{p, /, c} m+1 )C w ({p, /, c} m+1 ) 

J (m + 1)! 

X (fa/}m+l|-R)fa(fa/}m+l) 


= ofaF] + d m+ i[F] , 


(A.4) 


where the action of the operator fi is given by 

(fa /falfaPm) = fafa /}rn)(fa , (A.5) 


with 


m N i i (fa/}m|pm) , (fa, /}m|I(*o) + K(t 0 ) + P|Pm) , A 

fafa/}m) = l + —-—7HT7 +--- , , ,ns - • (A.6) 


(fa/}m|Pm) (fa/}m|Pm) 

The density matrix for the real subtracted cross section, (fa /} m +i|fa, is defined as 

POO _ 

(fa/Wil-R) = (fa/}m+i|pS+i) - / dr (fafaWilfacWi) 

•J to r _ _ / _ . /1 


{s,s',c,c'} 


[d{p, /, s', c', s, c} m ] ^ , ( 0 w s, s (A.7) 

-i-(fa,/,« ,c,s,c} m |pW) T>(fa,/} m ) 


mi 


X (fa, /, fa fa s, c} m _|_i|7far)|fa, /, fa fa s, c} m ) , 


where we specified explicitly the full quantum correlations in the subtraction terms required 
for the removal of kinematic singularities. 7ii(t) is the Nagy-Soper real splitting operator 
without any approximations. The inclusive jet functions ifafa, f} m ) define the generation 
cuts. 

Let us now study the accuracy of Eq. (A.4) after expansion to next-to-leading order. 
The parton shower evolution operator, currently available in DEDUCTOR, reads 

U(t F ,t 0 ) = N LC+ (t F ,t 0 ) + [ dr U(t F , t )'H^ C+ (t)N lc+ (t , to) , (A.8) 

Jt 0 

where 'H I (r) denotes the spin averaged splitting operator in the LC+ approximation and 
N LC + (tp, to) the corresponding Sudakov form factor, as discussed in Section 2.6. Inserting 
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Eq. (A. 8) in a m [F] and expanding the evolution yields 


Gm[F] = —7 [ [d{pj ,c} m ]{F\{pJ ,c}ra)C w ({pJ ,c} m ){{p, f} m \n\p { F ) )F I {{p, f} m ) 

771. J 

+ —i f[d{p,f,c} m \ [ dr (F\T~L^ C+ (t)— V lc+ (t)\{p, f,c} m ) (A.9) 

ml J Jt 0 

x C w ({pJ,c} m )({pJ} m \p$)Fi({p,f} m ) + 0(a 2 s ) . 

Using the approximations 


V{t) = V lc+ (t) + 0(1/N 2 ) 


and 


C W ({P, f, c} m )({P, f }m|Prn' 1} ) = ({P, /> c}n*|Pm 1} ) + °(. 1 / N c) > 
Eq. (A.9) reduces to 


(A.10) 

(A-11) 


ov 


= [[d{pJ,c} m ](F\{pJ,c}m)({pJ,c}m\ \p i n}) + \p { ^) + F\p { F ) ) Pj({P>/}m) 

Tfl. J L J 

[d{p,f,c}m] [ dr (F\H L 1 C+ {T)\{pJ,c} m ){{p 1 f,c} m \p ( F > )F I {{p,f} m ) 

Jtn 


ml j Jt o 

+ C»(^,l/iV c 2 ,Acr) . 

(A.12) 

Here, Aa represents the error resulting from the limit tp —>• oo. 

Let us now turn to the cr m +i[F]. Expanding the shower evolution and making use of 
Eq. (A.11) to recover the leading colour correlated real matrix element, we find 


<Jm +1 [f J = 


[F| = / 


[ d { p , f , C } m+ j] 

(m + 1 )! 


{F\{p, /, c} m+ i)F/({p, /} m+ i) 


({P> /> c }m+l|pS+l) 


/to 


/ 


dr E ({s , C } m _)_i|{s, c} m -|_i) C^dP; /, c}m+l) 

{s,s',c',c}m+l 
[d{p,/, s',c',s,c} r 


m! 


■({p> /, s', c', S, c} m |Pm) E>({p, /}m) 

X ({p,/,s / ,c / ,s,c} m+ i|^/(r)|{p,/,s',c,s,c} m ) +0(a 2 ,l/A c 2 ) 


(A.13) 


= ^+i[^] - <+i[F] + 0{pii, 1 /IV 2 ) • 

In the following we will focus on the subtraction terms a^ +1 [F]. The summation over 
{c} m+ i can be eliminated by 


J[d{p, /, c} m+ i]C w ({p, /, c} m+ i)(E|{p, /, c} m+ 1 ) 

= f [d{p,f} m+1 ](F\{p,f} m+1 ) + 0(1/N 2 ) . 


(A.14) 
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By a further approximation to the colour correlator 


({c'WilMm+i) = ({c} m+ i|{c} m+ i) 5 ({c'} m+ 1 ; {c} m+ i) + 0{l/Nl) , (A.15) 


we recover leading colour correlations 

J[d{p,f}m+ 1] ( F \{P’f}™+l) ({c}m+l\{c}m+l) 

{5}m + l (A.16) 

= J[d{p,f,c} m +i](F\{p,f,c} m+ i) . 

In the following we rename c —>• c and drop the second colour index in (fa, /, s', c, s, c\ m +i\- 
In addition, we use 

ftj(r) = Wf c+ (r) + 0(1/N 2 ) , (A.17) 

to obtain the following form for cr^ +1 [_F] 

Cm+lI-F] = ( m + l)\ m \ j [4p. /i 4m+l][ rf {p. /. *'> c', S, c}m]Ff({p, /} m )F;({p, /} m+ i) 

poo _ 

X fa|fa,/,c} m+ i) / dr ^ ({s'} m+ i|{s} m+ i) ({p,/,s',c',s,c} m |/9^ ) ) 

to {s,s'} m+ ! 

x ({p,f,s',s,c} m+1 \nj C+ (r)\{p,f,s',c',s,c} m ) + 0(a 2 ,l/N 2 ) . 

(A.18) 


It is now necessary to remove spin correlations. This can only be achieved, if (F\{p, /, c} m +i) 
and Fj({p, /} m +i) are sufficiently inclusive to allow an azimuthal average. A typical case 
would be an observable with NLO accuracy, where the azimuthal average corresponds to 
the complete phase space integral over unresolved jets. On the example of top-quark pair 
production with at least one jet, this would be any observable insensitive to additional jets. 
Once an azimuthal average has been performed, spin correlations vanish as demonstrated 
in Ref. [16] (Chapter 12) 

[ X] (fa'jm+llfalm+l) {{pJ,s',S,c} m+1 \n^ C+ (T)\{p,f,s',c!,S,c} m ) 

{s,s'} m+ l (A. 19) 

= ({S , }m|{s}m) (fa, f,c}m + l\ / H^ C + (T)\{pJ,c',c} m ) , 


which can be further transformed with 

(•fa, /, c}m+l|^/ C+ fa)lfa, /, C', C} m ) = (fa, /, C} m+ i ° + (t) |fa, /, £} m ) 

x6({c'} m -{c} m )+0(l/N 2 ) . 

Substituting Eq. (A. 19) and Eq. (A.20) into Eq. (A. 18) we obtain 

< 4 + l[F] = + 1 1 ) !m! J{ d {P’f’ C }rn+l][d{pJ,c} m ]{F\{p,f,c} m+ i) 

x f dr({p, /, c} m+ i|^j C+ (r)|fa,/, c} m )F/(fa, /} m )F/(fa,/} m+ i) 

Jto 

X X (fa'}m|fa} m ) (fa, /, s', s, c} m \p$) + 0(a 2 , 1/N 2 ) . 

{•M'l™. 


(A.20) 


(A.21) 


48 - 




The last line of Eq. (A.21) corresponds to the leading colour, spin averaged quantum density 
matrix 

({p, f, s', S, c} m \p^>) = ({p, /, c} m \p^) . 

With these approximations the final expression for reads 

°m+l[ F ) = 7- An, I [ [d{pj,c}m+l][d{p,f,c} m ] f dr ( F\{pJ,c}m+l ) 

(m + l)\m\J J to 

X ^>({P, /}m+l)({P, f, c}m+l\'H I [ C + {t)\{p, f, c} m ) 
x {{pJ,c}m\p$)Fj({p, f} m ) + 0(a 2 s ,l/Nc,Aa) . 

Combining a m [F] from Eq. (A.12) with cr^ +1 [F] and <r^ +1 [F] from Eqs. (A.13) and (A.23) 
yields 


(A.22) 


(A.23) 


= er 


+ 


a[F\ = a m [F} + a* +1 [F} - a^ +1 [F] 
r LC—NLO 


1 


[d{p,f, c} m+1 ][d{p, /, c} 


Ff 

mi I dr (F\{p,f,c} m+1 ) 
to 


(m + l)!m! 

x ({p, /, c} m+ i|-Hj C + (r)|{p, /, c} m )({p, f, c} m \p^) 
x F/({p,/} m ) [l-E/({p,/} m+ i)l +0(a 2 s ,l/N 2 ,Aa) 


(A.24) 


The last term vanishes for suitable generation cuts as specified in Section 3.3. In conclusion, 
we have shown that the simplified matching approach yields a prescription that is accurate 
at NLO to the level of leading colour. 
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